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ABSTRACT 

We investigate in detail the white dwarf cooHng sequence of the globular cluster Messier 4. 
In particular we study the influence of various systematic uncertainties, both observational and 
theoretical, on the determination of the cluster age from the white dwarf cooling sequence. These 
include uncertainties in the distance to the cluster and the extinction along the line of sight, as well 
as the white dwarf mass, envelope and core compositions and the white dwarf-main sequence mass 
relation. We find that fitting to the full two-dimensional colour-magnitude diagram offers a more 
robust method for age determination than the traditional method of fitting the one-dimensional 
white dwarf luminosity function. After taking into account the various uncertainties, we find a 
best fit age of 12.1 Gyr, with a 95% lower limit of 10.3 Gyr. We also perform fits using two 
other sets of cooling models from the literature. The models of Chabrier et al (2000) yield an 
encouragingly similar result, although the models of Salaris et al (2000) do not provide as good 
a fit. Our results support our previous determination of a delay between the formation of the 
Galactic halo and the onset of star formation in the Galactic disk. 

Subject headings: globular clusters; individual (Messier 4), age - stars: white dwarfs, luminosity function. 
Population II - Galaxy: halo 

1. Introduction 
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The use of globular clusters as a laboratory to 
study stellar evolution has a long and venerable 
tradition. In particular, a comparison of theoret- 
ical models with the observed properties of clus- 
ter stars allows one to determine an age for the 
cluster and the use of these ages as a lower limit 
to the age of the universe has become a standard 
cosmological test. The age determination rests 
on fitting models to the observed location of the 
main-sequence turnoff in the cluster colour magni- 
tude diagram. Significant effort from many groups 
(reviewed in VandenBerg, Bolte & Stetson 1996; 
Krauss & Chaboyer 2003) has been invested in ex- 
amining the detailed dependence of this method 
on both observational and theoretical systemat- 
ics. The current state of the art is summarised 
in Krauss & Chaboyer (2003), and yields a 95% 
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range for the age of the Galactie cluster system of 
^^■^-2.2 Gyr. The principal observational limita- 
tion is the uncertainty in the determination of the 
distance to a given cluster (e.g. Bolte & Hogan 
1995) while there are also several theoretical un- 
certainties, which include the treatment of con- 
vection, the treatment of helium diffusion and the 
conversion from effective temperature to observed 
colour. 

A direct consequence of the large age of the 
globular clusters is that a significant fraction of 
the stars that were originally born in the clus- 
ter have completed their stellar evolution, ending 
their days as white dwarf stars. A detailed study 
of the white dwarf population of a cluster allows 
us the opportunity to probe the upper main se- 
quence of a population II system as never before, 
as well as providing the opportunity to determine 
a cluster age by an entirely different method. As 
with the main-sequence turnoff method, there are 
both observational and theoretical systematics to 
be considered, but they are different systematics, 
so that we may hope to learn independent informa- 
tion by comparing the results of the two methods. 

The problem with studying white dwarfs is that 
they are faint and difficult to detect in the crowded 
environment of a globular cluster. While the clus- 
ter main sequence has been studied since the 1950s 
(Arp, Baum & Sandage 1952; Sandage 1953), the 
detection of cluster white dwarfs had to await the 
advent of the Hubble Space Telescope (Elson et 
al 1995; Paresce, de Marchi & Romaniello 1995; 
Richer et al 1995). To obtain useful age informa- 
tion requires extending the observed cooling se- 
quence to very faint magnitudes, requiring a sig- 
nificant investment of HST time. In Richer et al 
(2002) & Hansen et al (2002, hereafter paper I) 
we reported the first investigation of a large num- 
ber of cluster white dwarfs in the nearby globular 
cluster Messier 4. Comparison of the observed lu- 
minosity function with a single predetermined set 
of cooling models yielded an age for the cluster 
12.7 ± 0.7 Gyr (2cr random error). 

This paper represents an extension of the anal- 
ysis performed in paper I. The purpose is twofold. 
As this is the first data set to reach these magni- 
tude levels in a population II system, we investi- 
gate some general features of the observed cooling 
sequence and the physics that underlies the dis- 
tribution of white dwarfs in the colour-magnitude 



diagram. The second goal is to investigate in detail 
the various systematic effects, both observational 
and theoretical, that influence the age we infer for 
the cluster from the cooling sequence. 

The organisation is as follows. Section 2 de- 
scribes our definition of the white dwarf cooling 
sequence and the derivation of physical proper- 
ties. Of particular importance is our treatment of 
the extinction along the line of sight. Section 3 
then describes the details of the theoretical mod- 
els which we use to fit the cooling sequence. In 
section 4 we discuss some general features of the 
cooling sequence, noting the location of important 
physical regimes, before we proceed to the detailed 
model fitting in Section 5. The reader impatient 
for the bottom line may proceed to Section 6. 

2. The White Dwarf Sequence 

The data are based on two HST programs, 
GO 5461 and GO 8679, taken in cycles 5 and 9 
respectively. The data consist of 15 x 2600s expo- 
sures in F555W and 9 x 800s exposures in F814W 
in the first epoch (Richer et al 1995, 1997; Ibata et 
al 1999) and 98 x 1300s in F606W and 148 x 1300s 
in F814W in the second epoch (Richer et al 2002, 
2004a; Hansen et al 2002). Here we review briefiy 
the data reduction procedure - more extensive de- 
scriptions can be found in Richer et al (2004ab). 

The data were preprocessed according to the 
procedure outlined in the ALLFRAME cookbook 
(Stetson 1994; Turner 1997). For each individ- 
ual frame, we used DAOFIND with a A-a find- 
ing threshold to generate a star list. The reg- 
istered frames were then combined into 16 im- 
ages, corresponding to the two filters in each of 
the two epochs for each of the four CCDs in the 
WFPC2 camera. The registered frames were com- 
bined using IRAF's IMCOMBINE task to reject 
the brightest n pixels, and averages (not medians) 
were constructed of the remaining pixels. This 
method of combining was preferred as it provided 
us with higher S/N images; the noise in the av- 
erage of two frames will be similar to the noise in 
the median of three frames (Stetson 1994). Hence, 
provided we reject less than 1/3 of the pixels, the 
noise in the high-pixel rejected averaged frame will 
be less than that in a median of all the data. 
In practice this technique allowed us to use 93% 
of the FQOQW frames in the long exposures in 
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our field and 94% of the FSIAW images in con- 
structing the final images. Finally, we generated a 
finding list for the combined frames in the second 

epoch data. 

The M4 fields are dominated by cluster stars. 
As a result, the transformations used to register 
the frames were dominated by cluster stars (in- 
deed, as the matching radius was decreased, the 
non-cluster stars were preferentially rejected). A 
consequence of this is that non cluster stars ap- 
pear to move when frames from the two epochs are 
blinked. To improve transformations between in- 
dividual frames, we generated a list of stars likely 
to be cluster members by isolating those near the 
origin on the {x„ew - Xom) vs. {y„ew - Void) di- 
agram (membership is confirmed by the narrow 
locus such stars trace out in a CMD). The list of 
candidate cluster stars is on the coordinate sys- 
tem of the transformed frames. We then per- 
formed the inverse transformation of this coordi- 
nate list back onto the system of each of the pre- 
registered frames. Using this new list the individ- 
ual (pre-transformed) frames were re-reduced us- 
ing the ALLSTAR program with the recenter op- 
tion enabled. These 'cluster only' ALLSTAR files, 
which have only cluster candidates, were matched 
with DAOMASTER, using a 20 term transforma- 
tion equation. By experimentation, we found that 
the 20 term transformation reduced systematic er- 
rors. Transformations with fewer terms lead to 
apparent streaming motions in plots of dx vrs dy, 
indicative of poor fits - often occurring near the 
edges of the chips. As prescribed by Stetson, the 
DAOMASTER matching radius was reduced un- 
til it was a few times larger than the standard 
deviation of the differences in the transformed po- 
sitions. The 20 term transformation between the 
frames (based on cluster stars only) was used to 
produce new transformed images. Note however 
that this time we used the M0NTAGE2 program 
to expand the frames by a factor of 2 (pc CCD) 
or 3 (wf CCDs). The expanded, registered images 
were then combined again into four images, one 
for each filter in each epoch. 

By expanding the frames, we found that the 
measured radial dispersion in the cluster proper 
motion was reduced from 3.9 to 2.5 mas/year. The 
intrinsic dispersion amongst the M4 stars is esti- 
mated to be about 0.5 mas/year, so the measured 
dispersion is representative of the errors associ- 



ated with our technique. In effect, by expanding 
the frames we implemented a form of 'drizzling' 
(Pruchter & Hook 1997). 

2.1. Photometry 

Photometry on the expanded frames was per- 
formed using ALLSTAR. Stetson's library PSFs 
were not appropriate for the expanded frames, so 
new PSFs for these frames were built following 
standard techniques outlined in the DAOPHOT 
manual. All objects in the master coordinate list 
were visually classified (as to stellar or extended, 
isolated or possibly contaminated). 

Photometry requires some care in this field be- 
cause the crowded field results in a high back- 
ground. To define the photometric zero-point, we 
define a set of fiducial stars in the least crowded 
parts of the field and perform aperture photome- 
try and use the synthetic zero-points from Holtz- 
mann et al (1995) to place them on the WFPC2 
system. All other stars are measured relative to 
this fiducial set of stars. 

To perform a proper comparison, we further- 
more need to convert the WFPC2 magnitudes into 
the ground system, using the transformations of 
Holtzmann et al. This is necessary as we hope to 
later compare to various other models in the lit- 
erature, most of whom are only available in the 
ground-based system. The transformations are in 
terms of V-I, whereas we really want them defined 
in terms of F606-F814. To that end, we determine 
the reverse transformations 

V-I = -0.018 + 1.617 {F606W - FSIAW) 
-0.066 {F606W - FSUWf (1) 

and 

V = F606VK + 0.014 + 0.344 (^'GOOVF - ^8141^) 

-^0.037 {F606W - FSUWf (2) 

for F606M^ - F8UW < 0.66 and 

V = F606W - 0.018 + 0.407 {F606W - ^^814^^) 

-h0.015 {F606W - FSUWf (3) 

for F606W - F8UW > 0.66. 

These transformations are important as our 
deep exposures were in the F606W bandpass, 
which is broader than the more usual F555W 



3 



bandpass (and therefore attractive for such expen- 
sive projects as this one) but has the consequence 
that there is a significant colour-dependence in the 
WFPC2/ground transformation. Note this also 
has an important interaction with the extinction, 
as the colour transformations are determined us- 
ing bright standard stars, whose spectral energy 
distributions may imperfectly mimic that of a red- 
dened white dwarf. We will return to this point in 
§ 2.4. 

2.2. Proper motions 

At this stage, the proper motion of a star (rel- 
ative to the cluster) can be determined by exam- 
ining its shift in x and y. Cluster stars can be 
isolated in the data as they will have small proper 
motions. An important point to note here is that 
we were matching long exposures in F606H^ (127.4 
ksec) and F814H^ (192.4 ksec) with much shorter 
exposures in (31.5 ksec) and F814W (7.2 

ksec) . The question thus arises as to how we could 
measure positions in the earlier epoch data for the 
faintest stars where the S/N is obviously poor. For 
a star at V = 29.0 the expected S/N (from the 
HST ETC) in our deep F606W frames is 6.7 while 
a star at this magnitude on the shorter F555W 
frames has S/N = 2.5. This is just about what we 
measured on our frames. The reason we were able 
to measure positions of such faint objects in the 
first epoch are severalfold. 

First, we applied the finding list from the deep 
frames to the shallower frames. If the background 
noise in an image is Gaussian, then a 2.5c7 posi- 
tive deviation will occur in about 6 pixels out of 
every 1000. If, in a 750 x 750 image (neglect- 
ing the vignetted areas at the low-x and low-y 
sides of the WFPC2 CCDs), we were to mark 
all of the 2.5(7 peaks as detected astronomical ob- 
jects, we would expect more than 3,000 false de- 
tections. However, if we consider only the area 
within 0.5 pixels radius of an object confidently 
detected in the long second-epoch exposure, we 
expect a probability ~ tt x 0.5^ x « 0.005 of 
finding a positive 2.5(7 deviation which is purely 
the result of random noise in the first-epoch im- 
age. That is to say, we expect of order 5 false cross- 
identifications for every 1000 correct re-dctections. 
Even at 1.8a, presumed re-identifications will be 
correct 19 times out of 20. For this reason, the 
knowledge that an actual astronomical object is 



present somewhere nearby based upon the long- 
exposure second-epoch images allows us to be con- 
fident that most of the claimed re-detections on 
the short-exposure first epoch images correspond 
to true re-detections. The first-epoch astrometric 
positions, then, while poorer than those of the sec- 
ond epoch, arc nevertheless good enough to distin- 
guish stars that are moving with the cluster from 
stars and galaxies that are not. 

Secondly, the proper motions of the M4 stars 
are actually quite large (about 1 HST pixel with 
respect to an extragalactic background over the 6 
year time baseline (Kalirai et al. 2004)). While 
the lower S/N image would not produce good 
enough photometry, the S/N is sufficient to give 
the centroid, which is crucial for astrometry. In 
Figure 1 we illustrate the quality of the proper 
motion separation between cluster and field from 
the deep images in the outer field. Some of the 
brightest stars do not exhibit clean separation due 
to their near saturation. However, it is clear from 
this diagram that most field objects can be eas- 
ily eliminated down to quite faint magnitudes by 
making the generous proper motion cut within a 
total value of 0.5 pixel (~ 8 mas/yr) of that of the 
mean cluster motion over the 6 year baseline of 
the observations. 

Isolating the cluster main sequence in this man- 
ner is easy, given the large cluster proper motion. 
To isolate the white dwarfs, a little more caution 
is required, as clearly our ability to centroid is ex- 
pected to get worse as the stars get fainter and 
so the proper motion errors will be larger for the 
faintest stars (which are of most interest to us). 
To determine the appropriate proper motion cut 
for our purposes, we will perform the separation 
in a different way than is usually done. First we 
apply a CTit in the colour-magnitude diagram, to 
isolate all the stars (both cluster and background) 
that are of appropriate colour and magnitude to 
be candidate white dwarfs. We then examine the 
proper motion displacements of only these candi- 
dates, in order to determine the appropriate sepa- 
ration into cluster and background objects, with- 
out being biased by the excellent separation of the 
brighter main sequence stars. 

Figure 2 shows the colour-magnitude diagram 
for all the stars in our sample, along with the 
colour cut we apply. Figure 3 shows the proper 
motion displacement (relative to the bright clus- 
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ter main sequence stars) of all the stars that lie 
blueward of this criterion. Clearly we can distin- 
guish cluster from background to very faint mag- 
nitudes. However, since we hope to obtain age 
discrimination by examining the structure in the 
luminosity function, there is a concern that the 
exact location of the cut could influence the finer 
details of what we call our observed luminosity 
function. Including a self-consistent treatment of 
the proper motion cuts in the incompleteness es- 
timates (§ 2.3) should take this into account but 
we will also further investigate this by imposing 
two different proper motion cuts to see if there is 
any significant difference. The cut designated A 
will be the default value and we will examine the 
influence of cut B on the results later. 

Figure 4 shows the resulting white dwarf cool- 
ing sequence for proper motion cut A. The total 
number of white dwarfs detected in the three WF 
chips is 272 (103 in WF2, 78 in WF3 and 91 in 
WF4). This is the data set we attempt to model 
in § 3. Note that henceforth we will operate in V 
and I, rather than F606 and F814. 

De Marchi et al (2003) have published an inde- 
pendent analysis of this data which reaches signif- 
icantly different conclusions from those presented 
here. In particular their luminosity function does 
not extend as deep as ours, causing them to con- 
clude that all one can assert is that the age of M4 
is larger than 9 Gyr. We consider their analysis to 
be flawed and disagree strongly with their result 
on several points, which we discuss in Appendix C 
and Richer et al (2GG4b). 

2.3. Incompleteness 

An important part of our modelling procedure 
in later sections is the transformation between the- 
oretical models and observational parameters. In 
a cluster environment, our ability to accurately 
measure the magnitude and colour of faint objects 
can be strongly affected by the presence of many 
brighter stars in the field. To accurately model 
the resulting errors we have performed extensive 
artificial star tests, adding stars of known colour 
and magnitude at specific locations in the field and 
then rereducing the frames to measure the result- 
ing dispersion in colour, magnitude and position 
of the artificial stars. 

The added stars were drawn from a fiducial 



model white dwarf sequence and 4 trials of 1250 
stars each were added into each chips in the field. 
The stars were added onto a grid which ensured 
that no artificial stars overlapped with each other. 
The position of the added star at the intersection 
points of this grid was given a small random off- 
set to insure that stars were not all added at the 
same positions in different trials, but the shift was 
small compared with the grid and PSF size. This 
technique (first used by Piotto & Zoccali 1999) 
allows large numbers of stars to be added into a 
single frame while ensuring that they do not in- 
terfere with each other (i.e. do not change the 
crowding statistics in the images). This guaran- 
tees that the corrections thus derived provide ac- 
curate incompleteness statistics. We kept track of 
all the stars, irrespective at which magnitude they 
were recovered. We also determined the position 
at which the star is detected - necessary to prop- 
erly model the effects of proper motion selection 
when the proper motion accuracy degrades with 
magnitude. From this we derived a "correction 
matrix" with input magnitude along one axis and 
recovered magnitude along the other, as a func- 
tion of proper motion cutoff. The diagonal ele- 
ments of this matrix are the recovery fractions for 
stars recovered within the same magnitude bin as 
input while the off-diagonal ones are the recovery 
fractions for stars recovered at magnitudes either 
brighter or fainter than the input value. Figure 5 
shows the total recovery fractions (within a proper 
motion cutoff of 0.5 pixels) for the three WF chips. 
Of course, this is only the crudest measure of the 
selection effects. For the purposes of modelling, we 
also required the distribution of recovered magni- 
tudes (and positions) as a function of input magni- 
tude. Figure 6 shows the Icr (i.e. 68% limits - the 
distributions are not strictly Gaussian) bounds of 
the distribution (for both stars recovered at mag- 
nitudes both brighter and fainter) as a function of 
input magnitude for the WF4 chip (the distribu- 
tions for the other chips are similar). Note that 
this distribution reflects only those artiflcial stars 
which were recovered within 0.5 pixels of their in- 
put position, so this should represent the scatter 
in the stars we are actually trying to model. 

2.4. Distance and reddening 

The principal uncertainty in the age determi- 
nations by main sequence fitting is the uncertain 
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knowledge of the distance and extinction, which 
influence directly the estimate of the turnoff mass 
and thus the age. Similarly, the uncertainty in the 
distance and extinction translates into an uncer- 
tainty in the masses of the white dwarfs, although 
the relation to the age is not quite as straightfor- 
ward, as we shall see. 

The distance and extinction quoted in papers I 
and II were determined by fitting low metallicity 
subdwarfs to the observed lower main sequence 
(Richer et al 1997). The quoted result was a 
distance modulus in V of 12.51 ± 0.09, with E(V- 
I)=0.51 ± 0.02. The reddening parameter for this 
line of sight is estimated to be a nonstandard 
Rv = 3.8, which yields a true distance modulus 
11.18, in good agreement with astrometric de- 
terminations (Peterson et al 1995) and Baade- 
Wesselink measurements (Liu & Janes 1990). 
These methods constrain the distance to M4 to 
be 1.73 ±0.14 kpc. 

While appropriate for most applications, a dis- 
tance and extinction determination using subd- 

warf fitting is a potential source of systematic er- 
rors in our case. For one thing, it depends on the 
assumed metallicity (and one of the attractions of 
the white dwarf age method is that it is suppos- 
edly not particularly sensitive to the metallicity). 
Furthermore, in practice one has to slightly adjust 
the colours of some of the local subdwarfs whose 
metallicities are close but not exactly the same as 
that of the cluster stars, thereby potentially intro- 
ducing errors from model atmospheres. Finally, 
there are some differences between the spectral 
energy distributions of main sequence stars and 
white dwarfs, so the final extinction in a given 
bandpass can be slightly different for the two dif- 
ferent classes. Errors in the extinction contribute 
to errors on the age determination, since the red- 
dening vector has a significant component parallel 
to the cooling sequence (Figure 4). For these rea- 
sons, we would like to determine the extinction 
directly for the white dwarfs. We do this as fol- 
lows. 

The apparent magnitude of white dwarfs (of 
fixed temperature) at the top of the cooling se- 
quence is determined by the distance, the extinc- 
tion and the white dwarf radius (and hence mass) . 
In fitting a model to the observed cooling se- 
quence, one can then trade off variations in these 
parameters while keeping the apparent magnitude 



fixed. As a consequence, different extinctions im- 
ply different masses at the top of the cooling se- 
quence if we insist that the model sequence overlap 
the observations. The mass dependence of white 
dwarf cooling means that the choice of extinction 
(and hence mass) has implications for the luminos- 
ity function at fainter magnitudes and thus affects 
the age determination. 

Thus we wish to explore the sensitivity of the 
relation between the inferred white dwarf mass 
and the extinction. At the hot end, the theoretical 
cooling and atmosphere models yield 

Mj,o = 2.52(l/-7)o + 11.48-51og (j^) 

where we have used the subscript '0' to denote 
unreddened quantities and R is the radius of the 
white dwarf, which we shall henceforth denote Rg 
(i.e. in terms of lO^cm) to distinguish it from the 
reddening parameter R. This fit is performed in 
the range 0.1 < {V — I)o < 0.4, corresponding to 
effective temperatures 7500 < Tes < 10'' K. This 
was chosen to be cool enough to avoid the depen- 
dence of the radius on the surface hydrogen layer 
mass (not negligible for the hottest white dwarfs) 
and hot enough to avoid the non-monotonic effects 
of molecular hydrogen absorption on the colours. 

We also determine an empirical fit to the ob- 
servations (transforming from the flight system 
to ground-based system via the prescriptions of 
Holtzman et al 1995) 

/ = 2.53(1/ - J) -h 22.27 (5) 

in the range 0.5 < ^ - / < 0.9. This is not the 
true top of the cooling sequence but rather the lo- 
cation of the observed cooling sequence that cor- 
responds to the approximate location of the tem- 
perature range defined by the preceding theoreti- 
cal fit. Note also that there is some subtlety in- 
volved in defining what we term the 'cooling se- 
quence', as we refer here to those white dwarfs 
that evolve from single star evolution. A hand- 
ful of outliers may lie slightly above this sequence 
due to the presence of binaries - either as a result 
of the truncation of stellar evolution in binaries, 
which leads to undermassive white dwarfs or sim- 
ply due to two unresolved white dwarfs of similar 
age. Such undermassive white dwarfs have been 
observed in the field (Bergeron, Safii'er & Liebert 
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1992), in clusters (Cool ct al. 1998) and indeed in 
this cluster (Sigurdsson et al 2003). We have thus 
excluded a handful of objects that lie above the 
readily apparent cooling sequence. An example of 
such a fit is shown in Figure 7. 

Now, I = Mi + ^ + Ai andV -I = {V - I)o + 
Ay — Aj. Setting A/ = a Ay, we can eliminate 
V — I and / to get 

Av = 0.96 - 0.025(y - /) - 12.26 log Eg (6) 

where we have used a = 0.6 (we have determined 

that this number is robust for 3.1 < R < 3.8 
by using model atmospheres, extinction curves 
from Fitzpatrick(1999) and bandpass models from 
Holtzman et al (1995)). If we perform the fit in the 
middle of the fitting range {V — I — 0.6), we have 
a relation between the extinction and the white 
dwarf radius (and hence mass) which must be sat- 
isfied to satisfy the observed location of the cooling 
sequence 

Av = 0.94 ± 0.1 -12.261ogi?9 + 2.45(5^ (7) 

where we have also added the effect of a distance 
uncertainty S/j,. Thus, we observe a strong covari- 
ance between Ay and Rx). Note also the weak de- 
pendence on colour, which supports our assump- 
tion that the white dwarf mass is approximately 
constant over this colour range. This method is 
similar to the method used by Renzini et al (1996) 
and Zoccali et al (2001) to measure the distances 
to the clusters NGC 6752 and 47 Tuc by com- 
paring the cluster cooling sequence with a set of 
nearby white dwarf standards. However, in our 
case the sequences arc shifted parallel to the red- 
dening vector, rather than vertically in the CMD. 

What is a reasonable mass? The models of Ren- 
zini & Fusi-Pecci (1988) predict 0.53 ± O.O2Af0, 
which Renzini et al (1996) argue is a robust pre- 
diction. Cool, Piotto & King (1996) find a mass of 
0.55 ± 0.05Mq for the top of the cooling sequence 
in NGC 6397, assuming a standard distance and 
extinction. Similar numbers have been derived for 
the standard values in M4 as well (Richer et al 
1997). These numbers are slightly smaller than 
the mean value in the solar neighbourhood (e.g. 
Bergeron, Saffer & Liebert 1992), where the dis- 
tribution is sharply peaked near M ~ O.6M0. 
This would require Ay = 1.7 for our nominal dis- 
tance, although an extinction Ay = 1.3 can be 



obtained with fi ~ 11.03, which is just within 
the observational errors. We will return to the 
issue of distance/extinction in § 5.10. In princi- 
ple, the brightest white dwarfs in globular clus- 
ters are now accessible to spectroscopic analysis 
with large telescopes from the ground (Moehler et 
al 2000; Moehler 2002), offering the prospect of 
direct spectroscopic gravity (and hence mass) de- 
terminations. Initial results for the hottest white 
dwarf in NGC 6397 find a rather low gravity 
{logg ~ 7.3) suggesting an undermassive white 
dwarf (~ O.36M0), characteristic of binary evo- 
lution (e.g. Kippenhahn, Kohl & Weigart 1967). 
This illustrates that a sample of at least several 
gravity determinations per cluster are necessary to 
properly isolate the characteristic mass emerging 
from single star evolution. At present the spectra 
for M4 are not good enough to measure gravities, 
but do show that the hottest dwarfs have hydrogen 
atmospheres (Moehler 2002). 

Thus, our procedure is as follows. For each 
given model, we choose a characteristic white 
dwarf mass at the top of the cooling sequence. 
This then determines the appropriate extinction 
for that particular model, given an assumed dis- 
tance. In this manner we ensure that each of our 
models provides a consistent fit to the top of the 
cooling sequence. However, for the purposes of 
defining a default model, we start with a value of 
0.55 Mq, which implies Ay = 1.39 for our nom- 
inal distance /x = 11.18. Although this value is 
very similar to that determined by the subdwarf 
method, it has been determined directly from the 
white dwarf sequence, so that the agreement re- 
flects an encouraging concordance. 

3. Theoretical Inputs 

We have discussed above the characterisation 
of the observed Messier 4 white dwarf cooling se- 
quence. In order to determine an age we must 
fit a theoretical model to the observed distribu- 
tion of colours and magnitudes. We have already 
described two essential elements of the model - 
namely the model of the observational errors and 
the determination of the extinction and reddening. 
In addition we now describe the various elements 
of the cooling models. 
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3.1. Atmosphere Models 

Atmospheric models arc nciccssaiy to convert ef- 
fective temperatures and luminosities into observ- 
able quantities such as colours and magnitudes. 
This is a major issue in the determination of ages 
using main sequence stars, as the influence of met- 
als on the spectrum can lead to colour variations 
depending on the exact metallicity of the cluster 
stars. White dwarfs, on the other hand, have pri- 
marily hydrogen and helium atmospheres and so, 
in principle, are less susceptible to this kind of sys- 
tematic error. Even the few white dwarfs which 
do show metals in the atmosphere show only trace 
amounts which are not enough to significantly in- 
fluence broadband colours. 

There is one atmospheric subtlety which needs 
to be considered in white dwarf atmospheres. For 
white dwarfs at sufficiently cool effective temper- 
atures (Tcff < 5000 K), atmospheric hydrogen has 
a significant molecular component - which results 
in coUisionally induced absorption by H2 which 
can have a significant influence on the colours 
(Mould & Liebert f 978; Bergeron, Saumon & We- 
semael 1995; Hansen 1998; Saumon & Jacobsen 
1999). We will adopt as our default Teflp-colour 
transformations those of Bergeron, Wesemael & 
Beauchamp (1995) (augmented where necessary 
by fits to the models of Bergeron, Leggett & Ruiz 
2001). We also use a model we calculated ourselves 
(Hansen 1998; 1999) but the Bergeron models pro- 
vide a better fit (§ 5.7). 

Figure 8 shows the fit of an atmospheric model 
to the observed cooling sequence, assuming an 
extinction Ay — 1.39 and a white dwarf radius 
R = 8.9 X lO^cm. The model is clearly an excel- 
lent fit to the main locus of points. The deviation 
at the upper end is because, at the very highest 
temperatures, there is a small dependence of ra- 
dius on temperature which leads to slightly lower 
magnitudes. The exact shape of a true cooling se- 
quence will depend on how the mass varies along 
the cooling sequence. 

3.2. Cooling Models 

In addition to the atmosphere models, which 
relate T^s to the observed magnitudes, we also re- 
quire a set of cooling models, which determine the 
rate at which each white dwarf cools. White dwarf 
cooling is, in principle, quite simple to calculate. 



The luminosity is powered by residual thermal en- 
ergy left over from the prior nuclear burning his- 
tory. It is even insensitive to the exact temper- 
ature at which the white dwarf forms (since hot 
white dwarfs cool rapidly by neutrino emission - 
which means that the slow cooling begins where 
neutrino cooling ends, for all white dwarfs). Al- 
though cool white dwarfs do possess convection 
zones on the surface, such zones contain negligible 
mass and the parameterisation of convective effi- 
ciency has little effect on the cooling rate (unlike 
for main sequence stars). 

However, white dwarf cooling models do have 
their own set of uncertainties. 

• White dwarfs possess a carbon/oxygen core 
surrounded by thin shells of helium and hy- 
drogen. The thickness of these surface layers 
has an effect on the cooling as they regulate 
the rate at which heat is lost from the star. 
Stellar evolution models provide some guide- 
lines as to the expected size of these lay- 
ers, and asteroseismology offers the promise 
of one day determining the layer masses di- 
rectly, but for now the mass of the hydro- 
gen and helium layers remain a free param- 
eter which we will vary within the expected 
range. Figure 9 shows the influence of differ- 
ent layer masses on the cooling ages for our 
default model defined below. 

• In addition to the surface layers, the inter- 
nal chemical compositional profile in the car- 
bon/oxygen core can influence the cooling 
by virtue of its influence on the global heat 
capacity. Clearly a white dwarf that can 
store more heat will cool more slowly. This 
C/0 compositional profile is a function of 
mass and is also affected by the behaviour 
of the main sequence models and by the 
^^C(a, 7)^^0 reaction cross section, which is 
poorly known. 

• Cool white dwarfs crystallise in their cores, 
resulting in a release of latent heat which 
provides an extra source of heat to slow the 
cooling. In addition, there is the possibil- 
ity that the equilibrium mixture of carbon 
and oxygen may be different in the solid and 
liquid phases. In that case the resulting ad- 
justment in the density structure can lead to 
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an additional release of energy termed chem- 
ical separation energy. This too will depend 
on the nature of the compositional profile, as 
well as the phase diagram of the C/0 mix- 
ture. 

In subsequent sections we will investigate the 
influence of these various effects, both by including 
variations in our own set of models and by making 
comparisons with others in the literature. 

3.3. The Main sequence lifetime 

In order to get a proper estimate of the true age, 
we need to associate a given white dwarf with a 
main sequence progenitor, as the true age of the 
star is the sum of the white dwarf cooling age and 
the additional main sequence lifetime associated 
with that object. The main sequence lifetime is 
affected by overall metallicity and so we must be 
careful to choose an appropriate model and not 
simply appropriate such fits as have been used in 
previous papers on determining ages from white 
dwarf cooling models. We will use the li&^time- 
mass relation from the models of Hurley et al 
(2000) for Z = O.IZq. Figure 10 shows the re- 
lation between mass and pre- white dwarf lifetime 
for models of solar metallicity, and models of lower 
metallicity. 

3.4. The Main Sequence - White Dwarf 
Mass relation 

A related quantity is the relationship between 
main sequence mass and white dwarf mass. In 
principle, this can be predicted from a set of stel- 
lar models as well. However, our incomplete un- 
derstanding of mass loss in the advanced stages 
of stellar evolution makes this a less than secure 
proposition, and we also make use of empirically 
determined prescriptions. This is traditionally 
done by measuring the mass of white dwarfs and 
main-sequence turnoff stars in a variety of clusters 
of different age (Weidemann 1987, 2000; Koester 
& Reimers 1996; Claver et al 2001). Figure 11 
shows the range of relations we will explore in this 
paper, compared to the empirical relation of Wei- 
demann (2000). The solid curve denotes our de- 
fault relation, based on that of Wood (1992) 

M^d = O.55M0exp (0.095(M„, - Mto)) (8) 



where AIto is the turnoff mass appropriate to a 
particular age. This has the same form as the 
relation A used by Wood, but we have adjusted 
it so that the turnoff mass always gives 0.55Mq 
for all ages. This is necessary to be consistent 
with our method of determining the extinction in 
§ 2.4. Other theoretical relations are those due to 
Dominguez et al (2000) and Hurley et al (2000). 
Some scatter is evident for main sequence masses 
> 2AIq. This is also approximately the amount of 
scatter accomodated by the observations (Claver 
et al 2001). We should note that most empiri- 
cal measures arc for solar or super-solar metallic- 
ity stars and that some difference is expected for 
lower metallicity stars. We will test this by testing 
predictions for both solar and sub-solar metallicity 
models. 

One encouraging point to note is that the scat- 
ter in the predictions is primarily for main se- 
quence masses > 2M0, which means that cluster 
ages > 10 Gyr are not significantly affected. We 
will return to this point in § 5.4. 

4. Features of the cooling sequence 

Given that this dataset is the most extensive 
sample of population II stellar remnants known, 
it is of interest to discuss a few general features 
of the cooling sequence before turning to detailed 
model fits. Along our cooling sequence we can 
identify several important epochs in the life of a 
white dwarf. 

4.1. Hot white dwarfs 

Young, hot white dwarfs cool very rapidly - 
the hottest cool by neutrino emission. The neu- 
trino contribution to the cooling rate drops to less 
than 10% within 3 x 10^ years. The effective 
temperature at this point is ^ 25,000 K, which 
corresponds to approximately V — I ^ 0.29 and 
/ ~ 22.9 for this cluster. Our cooling sequence 
has 3 stars that lie above this point, suggesting 
the current birthrate of white dwarfs in our clus- 
ter field is ~ 10"'' yr~^. Rather than rely on just 
three stars, we can also estimate this rate by deter- 
mining the cooling age of the tenth faintest white 
dwarf on our sequence. This occurs at / = 24 and 
corresponds to T^s ~ 9700 K for our default pa- 
rameters. The cooling age for a 0.55 M© dwarf 
at this temperature is 5 x 10^ years, yielding a 
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rate estimate ~ 10/5 x 10^ 2 x IQ-'^yr'^ If 
the current white dwarf birthrate is this high, we 
expect ~ 10^ X 2 X 10"'' ~ 20 Horizontal Branch 
stars in our field. Richer ct al (1997) estimate 
only ~ 2, based on scaling the relative numbers of 
horizontal branch and giant stars. However, the 
cluster shows evidence for significant mass segre- 
gation. We have fit multi-mass Michie-King mod- 
els to the full M4 data set (Richer et al 2004a), 
including data from fields at smaller radii. These 
fields are also proper motion selected, using first 
epoch data from the program GO-5461 and second 
epoch archival data from program GO-8679. On 
the basis of these models we can relate the mass 
fraction in our field to the global value and can 
thus convert the birthrate in our field to a global 
white dwarf birthrate ~ 1.5x10"^ year. This now 
agrees well with the expected birthrate given the 
total number of observed Horizontal Branch stars 
in the cluster (~ 150/lOVs ~ 1-5 x IQ-^yr'^). 
We will return to the issue of mass segregation in 
§ 5.2. 

4.2. DA versus DB white dwarfs 

In the Galactic disk, some fraction of white 
dwarfs show helium atmospheres (DB) versus hy- 
drogen atmospheres (DA). We will discuss this 
in the context of the cool models later, but for 
the purposes of the upper part of the cooling se- 
quence, it is important to note that the choice of 
atmospheric composition can result in a change 
^ 0.2 magnitudes at fixed colour. This can ex- 
plain much of the dispersion observed in the up- 
per part of the cooling sequence but is too small to 
account for the 'overluminous' white dwarfs which 
lie above the main sequence. We will identify 
these as potential binary members below. For hot- 
ter white dwarfs, hydrogen atmospheres are much 
more common than helium atmospheres (e.g. Sion 
1984), so we use hydrogen atmosphere models to 
determine the extinction in § 2.4. This is further 
justified by the detection of hydrogen in the first 
spectra of hot white dwarfs in M4 (Moehler 2002). 

4.3. The ZZ Ceti strip 

As the white dwarf cools, the effective tem- 
perature begins to drop below the point where 
the atmospheric constituents remain ionized. As 
the white dwarfs pass through the range Teff = 
11,500-12,500 K, hydrogen atmospheres acquire 



a thin convection zone near the surface. It is 
also in this range of effective temperatures that 
we find the ZZ Ceti stars. These are stars which 
display photometric variability due to the presence 
of pulsations. The nature of the extinction uncer- 
tainty means that the blue edge could lie in the 
colour range V — I = 0.35 0.53, and the red edge 
V — I =0.42-0.60. This range encompasses about 
7 white dwarfs on our CMD, although only one or 
two would actually be ZZ Ceti stars (depending 
on where it is truly located). A study of the vari- 
ability in this field (Ferdman et al 2003) revealed 
no statistically significant photometric variations 
among these stars, but this does not pose a strong 
constraint as the integration times (1300s) are of 
order the variability timescale for ZZ Ceti stars, so 
that much of the variation would be washed out. 

4.4. The Luminosity Function Jump 

The most striking feature of the cluster white 
dwarf luminosity function is the sharp jump in 
the numbers for 7 > 26. This is a direct conse- 
quence of the behaviour of the white dwarf cool- 
ing curve. As the white dwarf continues to cool 
beyond the ZZ Ceti strip, the opacity of the now 
neutral atmospheric hydrogen decreases. The pho- 
tospheric pressure P (x g/n and the resulting in- 
crease in the pressure affects the outer boundary 
condition of the cooling models, causing the cool- 
ing to slow down. This results in a pile-up of stars 
at this point - as shown in Figure 12 (for a O.6M0 
model). The approximate location of this jump 
is at T^id ~ 5 Gyr and logL/L© ~ —4.1. It is 
important to note that this feature is primarily a 
result of the effective temperature - thus it does 
not depend strongly on the mass or core composi- 
tion of the white dwarf. This in turn implies that 
it should be a universal feature of all models, so 
that there is little age information in the location 
of the jump. 

4.5. Core crystallisation 

Crystallisation of the core material occurs at 
roughly the same time as the change in the outer 
boundary condition. Thus, the release of latent 
heat can also contribute to a temporary slowing of 
the cooling. However, there is an important dif- 
ference between these two potential contributors. 
The point at which crystallisation begins depends 
on the central density and thus the mass of the 
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white dwarf - as a consequence, white dwarfs of 
different mass receive the release of latent heat at 
different epochs in their cooling curves. However, 
the change in the atmospheric boimdary condition 
is a direct consequence of the effective tempera- 
ture and thus occurs at the same colour (except 
for changes in the reddening) for all white dwarfs. 
Crystallisation is very important, however, as it 
makes the cooling mass dependent (Figure 13). 
This mass-dependent behaviour contributes signif- 
icantly to differences in shape between the vari- 
ous model luminosity functions and which leads 
to much of the age discrimination in the methods 
we describe in § 5. 

4.6. Colour Evolution and the Luminosity 

Function cutoff 

At the faint end of the cooling sequence there 
are a couple of features one might hope to see. 
One is the change in the white dwarf colours due 
to the absorption of infra-red light by molecular 
hydrogen collisionally induced absorption (Mould 
& Liebert 1978; Bergeron, Saumon & Wcscmael 
1995). This causes white dwarf colours (at least 
in the optical and near-IR) to get bluer as they 
cool, instead of redder (Hansen 1998; Saumon & 
Jacobsen 1999).^ As shown in Figure 8, we cannot 
claim evidence of this distinctive feature, although 
such a statement is dependent on the assumed ex- 
tinction value. We also do not claim to have mea- 
sured a real truncation in the white dwarf lumi- 
nosity function - our model fits make use of the 
full distribution in colour and luminosity at the 
faint end to get constraints on the age. This point 
will be discussed at greater length in § 5. 

4.7. White Dwarf Binaries 

The binary fraction amongst white dwarfs in 
clusters is of interest for several reasons - not least 
of which is the suggestion that they might be ef- 
ficient producers of type la supcrnovae (Shara & 
Hurley 2002). Binary white dwarfs are identified 
as being overluminous either as a result of the 
light of two unresolved white dwarfs being inter- 
preted as a single object or else because the white 

^This has been referred to by some as the 'blue hook', which 
could lead to confusion in the globular cluster context, as 
this term is already used to describe a phenomenon associ- 
ated with anomalous Horizontal Branch stars! 



dwarf is of lower than average mass (and hence 
larger radius) due to the truncation of stellar evo- 
lution in binaries (Kippenhahn, Kohl & Weigert 
1967). This is a phenomenon frequently encoun- 
tered in globular clusters (Hansen, Kalogera & 
Rasio 2003), including this one (Sigurdsson et al 
2003). Note that there is an important difference 
between these two cases. For binaries containing 
two normal white dwarfs, only those of compa- 
rable cooling age will appear overluminous our 
estimate will only be a lower limit to the overall 
binary frequency. For binaries containing an un- 
dermassive white dwarf, no such correction applies 
as all binaries of this type are readily identifiable. 

We can get a rough idea by measuring the over- 
luminous fraction of white dwarfs in our field. Re- 
stricting our sample to / < 25 (to avoid thorny 
issues of incompleteness and photometric scatter) 
we find 5 overluminous white dwarfs out of a to- 
tal of 45 stars in Figure 4. There is also one 
anomalously massive white dwarf - which lies sig- 
nificantly below and to the left of the cooling se- 
quence. The identification of any single star as a 
cluster member is always open to question since 
some background objects can have fortuitous mo- 
tions which fall within the proper motion cut. 
However, this star has a proper motion well within 
the cluster cut and all measures are consistent with 
a cluster white dwarf. If it is a true cluster mem- 
ber then the location ~ 1.4 magnitudes below the 
cooling sequence implies a mass > IMq and an 
age ~ 10® years. If this object is the result of 
a merger of two O.5M0 dwarfs, then it suggests 
a minimum rate of white dwarf mergers ~ 10~® 
per year in this field alone. To convert this into a 
global rate we can again make use of the Michie- 
King models of Richer et al (2004a) , which suggest 
a global rate ~ 10 times larger than in our field 
alone, i.e. ~ 10"'^ yr~^.^^ 



'If the M4 rate is characteristic of Galactic globular clusters, 
then this implies a rate ~ I0^''yr^^ in the Galaxy as a 
whole. Such a number is of interest for the question of 
whether globular clusters are a significant source of Type la 
supcrnovae (e.g. Shara & Hurley 2002) . If this number is 
characteristic, then globular clusters contribute roughly 1% 
of the estimated Galactic rate. Of course, the rate could 
be higher in clusters more centrally condensed than M4 or 
indeed in the center of M4 itself (binary hardening due to 
external perturbations proceeds faster in the denser central 
regions). Also we note that we can only count white dwarf 
mergers that fail to become supernovae. However, given 
the characteristic distribution of white dwarf masses, our 
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The binary fraction of white dwarfs 0.11 ± 
0.05) we observe is interestingly similar to the one 
inferred for field white dwarfs (Holberg et al 2003) 
if wc restrict the sample to binaries that would re- 
main hard in M4 (The disk also contains a roughly 
similar fraction of wider white dwarf binaries that 
would be disrupted by three-body encounters in 
M4). This ratio is also consistent with our es- 
timates for the binary merger rate relative to the 
total birthrate of white dwarfs in our field. Our bi- 
nary fraction is slightly larger than the 1% inferred 
for the lower main sequence (Richer et al 2004a), 
but not necessarily a problem as the statistical er- 
rors are relatively large and the two numbers re- 
flect very different parts of the stellar population. 

Binaries are more massive than the average 
star, so our binary fraction estimate may not be 
appropriate for the cluster as a whole. The first 
epoch data (Richer et al 1995) for this cluster cov- 
ers a much wider area of the cluster than our sec- 
ond epoch data. Furthermore, a shallow second 
epoch taken by Bedin et al (2001) allows us to 
separate cluster white dwarfs by proper motion in 
the inner fields in the same manner as for our de- 
fault field, albeit to a much brighter magnitude 
limit. However, the depth achieved in these in- 
ner fields is sufficient to address the question of 
binary membership. Figure 14 shows the upper 
white dwarf cooling sequence in 3 fields closer to 
the core than our field. The data were divided 
into 3 radial bins so that roughly equal numbers 
of white dwarfs were found in each. We consider 
only the magnitude range 23 < F814 < 24.5, as 
for fainter white dwarfs the photometric scatter is 
too large and for hotter white dwarfs the temper- 
ature dependence of the white dwarf radius leads 
to an uncertainty. In all three fields, the num- 
bers are consistent with a 10% binary fraction, 
although this is based on small numbers. In the 
outer two fields we sec fractions 2/23 and 1/20. 
The inner field, however, shows a larger potential 
binary fraction of 0.31 ±0.15 (5/16). However, the 
photometric scatter in the inner field is larger, so 
that 3 of the binary candidates are uncertain, i.e. 
all fields are consistent with a 10% binary fraction. 

number is probably an upper limit. 



5. Modelling the Luminosity Function 

We have discussed above the various elements 
required to model the distribution of white dwarfs 
in the cluster. The model should match not only 
the distribution in luminosity but in colour as well. 
Since one can always keep a V or I band magnitude 
fixed while trading off mass/radius versus Toff, it 
is clear wc need to include colour information to 
get the best possible fit. Traditionally, white dwarf 
ages are determined by fitting models to the lumi- 
nosity function. In principle then, one could fit 
to both the V-band and I-band luminosity func- 
tions simultaneously. Alternatively, one could fit 
to the two dimensional distribution of stars in the 
CMD directly (this is referred to as fitting the Hess 
diagram in stellar population studies). The lumi- 
nosity function fits have the advantage that the 
one dimensional nature allows us to use smaller 
bins while still retaining enough stars in each bin 
to be statistically meaningful. However, only a 
two-dimensional fit can take proper self-consistent 
account of the colour information. 

To begin our model, we choose a (main se- 
quence) mass function of given slope x (where 

dN/dMms oc Mms'"^^^ ) in the white dwarf progen- 
itor mass range and sample it in a Monte Carlo 
fashion. Wc convert all main sequence stars to 
white dwarfs based on the initial-final mass rela- 
tion (our default model will assume relation 8) and 
assign a magnitude and colour based on the white 
dwarf cooling age once the main sequence age has 
been subtracted. This is then fed through the ap- 
propriate extinction and distance. Wc also Monte 
Carlo sample the results of our artificial star tests 
and assign an observed magnitude based on the 
probability distribution inferred for the observa- 
tional photometric scatter. 

• Fitting to the Luminosity F\inction: 

We then form V and I band luminosity func- 
tions from our model. We marginalise over 
the total normalisation to get the best fit 
(Avni 1976). The scaling is performed rel- 
ative to the number of stars in the ranges 
27.3 < V < 29.1 and 25.8 < / < 27.3. 
The fits to the observed luminosity func- 
tions arc applied over the entire magnitude 
range brighter than the 50% completeness 
limit (29.1 in V and 27.3 in I). The bins 
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arc 0.3 magnitudes large in both luminosity 
functions except for the top of the two cool- 
ing sequences, where the regions V < 27.3 
and / < 25.8 arc counted as one bin because 
of the low numbers of stars at these bright 
magnitudes. 

• Fitting the Hess diagram: 

A more quantitative method for taking ac- 
count of colour information is to bin the 
white dwarfs in the Hess diagram directly. 
Our grid for doing this is shown in Figure 15. 
For / < 26, we use a single bin, as the white 
dwarfs are relatively sparse and there is little 
age information in the fine structure of this 
part of the cooling sequence. After the lu- 
minosity function jump we use a finer grid, 
but use only those regions of the diagram 
that lie above the 50% completeness limits 
in both V and I. Once again, the overall nor- 
malisation is a free parameter over which we 
marginalise to get the best fit. Normalisa- 
tion is scaled relative to the total count in 
the bins shown in Figure 15. 

5.1. The Default Model 

We begin with a fit to our default model. The 
cooling models are based on the cooling models of 
Hansen (1999), with a hydrogen layer mass frac- 
tion qn = 10~^ and a helium layer mass frac- 
tion qHe = 10~^. We assume the white dwarfs 
at the top of the cooling sequence arc O.55M0 
and the initial-final mass relation as defined in 
equation (8). This choice of mass, along with the 
best fit distance (/io = 11.18) implies an extinc- 
tion Ay = 1.39. Colours and magnitudes are cal- 
culated using the transformations from Bergeron, 
Wesemacl & Bcauchamp (1996). The resulting fit 
is shown in Figure 16. 

The best fit to the Hess diagram is an excel- 
lent = 8.7 for 11-2=9 degrees of freedom (i.e. 
Xlof = 0.97), located at (a;, T)=(-0.85,12.4 Gyr). 
The 2a age Umit is T = 10.9-14.2 Gyr (where 
again, we interpret '2(t' as representing the range 
encompassed by the 95% confidence interval - the 
distribution is not gaussian). Strikingly, the lumi- 
nosity function fits yield somewhat different limits, 
with the V LF favouring older ages and the I LF 
favouring younger ages, although both fits overlap 
the 2(7 region of the Hess fit (albeit in different 



regions) . This is an illustration of why we need to 
examine the full colour information. An examina- 
tion of the 50% completeness limits in Figure 15 
shows that the V and I LF reflect different samples 
(since we exclude different stars on the basis of the 
completeness limits). Henceforth we will base our 
estimates on the Hess fits as this is the only way 
to treat the colour information self-consistently. 

At this point it is instructive to examine what 
distinguishes a good fit from a bad fit. Figure 17 
shows the Hess diagram for our best fit solution 
{x,T) =(-0.85,12.4 Gyr). Clearly both the colour 
and magnitude distributions are very well fit by 
the model. Figure 18 shows the V and I LF for 
the same solution. In both cases the shape of the 
luminosity function rise is well-modelled and the 
degradation in the is driven by the last few 
bins, which is just where we expect problems to 
occur when we try to model a two-dimensional 
distribution with a one-dimensional model. 

We also examine two cases of what we classify 
as a bad fit. Figure 19 shows the Hess fit if we 
retain the same mass function slope as our best 
fit model, but reduce the age to 10 Gyr. The 
Xdof = 5.7 for this model and is driven up by 
two problems. The first is that the model fails 
to reproduce the correct ratio of bright to faint 
white dwarfs, while the second problem is that the 
colour distribution of the model is too blue at the 
faint end. This is due to the fact that the lower 
age means the model white dwarfs are more mas- 
sive, hence fainter at late times, and so enter the 
faintest magnitude bins while still hotter and bluer 
than their less massive counterparts. Our second 
case of a bad fit is shown in Figure 20. In this case 
we examine what happens when we keep the same 
age as our best fit model, but increase the mass 
fimction slope x to x = 1.35 (Salpcter's value). 
In this case we again find the problem that the 
model cannot fit the ratio of bright to faint white 
dwarfs (although the problem is in the opposite 
direction in this case!) and furthermore, although 
the shapes of the colour distributions are roughly 
correct, the relative numbers of white dwarfs in 
the three different fainter magnitude ranges are 
no longer well fit. The Xdof = 2.5 in this case. 

The models used here are those that correspond 
most closely to that favoured by our prior prej- 
udices - the extinction is approximately that in- 
ferred from the subdwarf fits to the main sequence. 
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the white dwarf mass at the upper end corresponds 
to that expected from stellar evolution and the 
initial-final mass relation is similar to that used in 
other contexts. The cooling models also assume 
traditional hydrogen and helium layer thicknesses. 
We now wish to explore the sensitivity of our re- 
sults to a variety of model uncertainties. 

5.2. Number Counts 

Before we investigate other parameters, we note 
that there is a piece of information wc have not yet 
used. The mass function slope x that we quote is 
used only to determine the probability distribution 
of progenitor masses. At no point in the above 
analysis have we made any attempt to connect this 
to the observed number counts of white dwarfs 
and main sequence stars. In fact, wc will argue 
below that the relation between the quantity x and 
the true mass function is not trivial. However, to 
begin with wc will treat x as the true mass function 
slope. In that event, and assuming continuity of 
the cluster mass function, we can obtain additional 
information as follows. 

Each solution we obtain predicts a total num- 
ber of white dwarfs with F < 29. To proceed 
further, we need to make an assumption that con- 
nects the mass function of the white dwarf pro- 
genitors to that of the observed main sequence 
stars. The main sequence is successfully modelled 
(Richer et al. 2004a) as a power law mass func- 
tion with X ~ —0.9 over the range 0. 095-0. 56M0. 
To connect to the power law at the upper end, 
we assume continuity of dN/dM at some break 
mass M(,. Thus, integrating these two power laws, 
the ratio of white dwarfs to main sequence stars is 

Nms X2 * 0.56-^1 - 0.095-^1 ^ ' 

where X2 is the power law in the white dwarf pro- 
genitor mass range and xi is the corresponding 
lower main sequence value. The white dwarf pro- 
genitors range from the turnoff mass Mto to some 
mass Ml to be determined, which is the progen- 
itor mass of the faintest white dwarf we count. 
This is clearly going to be smaller than the canon- 
ical 8Mq limit for white dwarf production, as the 
white dwarfs that result from the most massive 
progenitors will have cooled beyond our detection 
limit and will not have been counted. We count 
520 main sequence stars in this limited range and 



222 white dwarfs brighter than V=29. However, 
the white dwarf numbers are too low because of 
incompleteness (the main sequence stars are some- 
what brighter and so not significantly affected by 
incompleteness). Once we correct for this, the ob- 
served numbers are equivalent to 290 white dwarfs 
brighter than V=29. Thus, with xi = —0.9, we 
find 

OM = -—M^'+°-^{M£'''-M^^') (10) 

X2 

Thus, given a choice of Mb, this yields an equa- 
tion for T in terms of x (since both Ml and Mto 
are functions of age T). This illustrates, in rough 
terms, the interrelation between x and T that we 
found in previous sections, although it is impor- 
tant to reiterate that the previous solution did 
not use any information about the lower main se- 
quence stars. 

There is still the matter of choosing M{,. Any 
mass Mb > 0.56Mo is allowed. Figure 21 shows 
the solution for three different choices of Mi, = 
0.6, 0.8, 1.OM0. They are compared to the (Hess 
fit) solutions generated for our default model from 
the previous section. The covariance of a; & T in 
the detailed solution is encouragingly similar to 
that obtained from this simple calculation with 
Mb ~ O.8M0. We emphasize that these values of 
X come from entirely different lines of argument. 

What is a reasonable range for xl Clearly a 
constant mass function slope x ^ —0.9 is not vi- 
able if extended into the regime of supernova pro- 
genitors or even massive white dwarfs. If a cluster 
loses a significant fraction of its mass during the 
course of stellar evolution, the cluster will dissolve 
(Chernoff & Weinberg 1990), so mass functions 
that are dominated by high mass stars arc not vi- 
able candidates. However, for ages of ^ 12 Gyr, 
Ml at our completeness limit is only ^ I.SMq. 
Clearly, there is little difference in postulating a 
break mass at O.SMq or at 1.2>Mq. In the latter 
case one can adopt the above solution without vi- 
olating the bounds of significant cluster mass loss. 

5.2.1. The Influence of Mass Segregation 

However, mass segregation is important in M4 
(Richer et al 2004a) and so the relation between 
the white dwarf number counts and the main se- 
quence number counts is no longer trivial. Michie- 
King models of the cluster suggest that the num- 
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bcr fraction of white dwarfs relative to all main 
sequence stars is approximately constant through- 
out the cluster. However, the relative numbers of 
white dwarfs to turnoff-mass main sequence stars 
(which more accurately reflect the spatial distri- 
bution of the white dwarfs when they had their 
progenitor masses) is a factor ^1.3 larger in our 
field than in the cluster as a whole. If we re- 
peat the calculation of equation (10) but extrap- 
olate from O.6M0 to O.8M0 (using .ti=-0.9) and 
count white dwarfs relative to only the mass range 
0.6 — 0.8Mq, then the number of main sequence 
stars in this range is 205 and thus N^d/Nms = 
1.414 in our field. This comparison is also a 
more robust measure of mass segregation because 
stars M > 0.6Mq are less susceptible to mass- 
dependent evaporation and tidal stripping. If we 
perform the calculation again using these numbers 
without correcting for mass segregation, we infer 
X2 ~ 0.4 for a representative case of T 12 Gyr 
(so that Ml ~ 1.3Mq and Mto = 0.87Mq) and 
Ml, = O.8Af0. However, if we divide N^d/Nms by 
1.3 to convert our field ratio to the global ratio, 
then we infer a value X2 ~ 1.4. 

This is only an illustrative calculation as it de- 
pends on several uncertain factors, but it shows 
that the effect of mass segregation is to make our 
value of X2 appear smaller than the true global 
value. The fact that a star loses mass in becom- 
ing a white dwarf means that we see more white 
dwarfs in our field than we would if they had re- 
tained the masses of their progenitors. Thus, we 
should be very cautious in our interpretation of 
X2 - it is primarily a parameter that describes the 
probability distribution of white dwarf progenitor 
masses in our field. Its relationship to the true 
global mass function is more uncertain and the 
true slope is likely to be somewhat steeper. In 
fact, King-Michie models for the full cluster sug- 
gest that the global value of a;2 ^ 1.3 (Richer et al 
2004a). 

5.3. Sensitivity to cooling models 

Using our newly defined default model, we can 
examine the systematic effects of the various un- 
certainties outlined in § 3. To keep track of the 
various effects, we record the results in Table 1. 



5.3.1. Hydrogen and Helium Layer Masses 

If we keep the helium layer mass fixed, and leave 
the internal core composition unchanged, but re- 
duce the surface hydrogen layer, we expect faster 
cooling. If we repeat the default model, now using 
a surface hydrogen layer mass qn = 10~^, then 
we obtain the contours of Figure 22. The best fit 
solution shifts slightly to (x, T) = (-1.6,12.4 Gyr). 
The lower age limit drops to 9.8 Gyr. However, 
this model has a higher than the default, so 
the lower age limit will be higher than this when 
we marginalise over the various model uncertain- 
ties. 

Similarly, there is an uncertainty in the he- 
lium layer mass. If we keep the default model 
the same (including qn = 10"**) but reduce the 
helium layer thickness to qne = 10"'^, we again 
expect the cooling to be more rapid. The solution 
becomes Figure 23. The shift in the 2a contours is 
small, although the best fit solution does move to 
(x,T) =(-0.5,11.7 Gyr). Again, the values go 
up across the board, indicating the best fit model 
at any given point is obtained with our default 
parameters. 

5.3.2. Internal Composition 

In addition to the uncertainty in the surface 
layer masses, there is also uncertainty in the com- 
position in the center of the star. The rela- 
tive amounts of carbon and oxygen vary through- 
out the star and are determined by the nuclear 
burning history of the star on the Horizontal 
Branch and AGB stages. The somewhat uncer- 
tain ^'^C{a,^y^O nuclear cross-section also plays 
an important part. Whereas the different sizes of 
the surface layers affect the rate of transport of 
heat through the star, the internal composition 
determines the overall thermal energy available to 
power the star, through both its influence on the 
heat capacity and extra sources of heat such as 
latent heat or chemical separation. 

The most conservative approach would be to 
simply compare with models of pure carbon and 
pure oxygen, but that is too naive, as no stel- 
lar evolution models predict pure compositions^^. 
Rather, we have calculated a second set of white 

^ After all, nobody calculates main sequence models with 
helium fraction=0! 
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dwarf cooling models using a different set of ini- 
tial conditions. Our default model was calculated 
using the chemical compositions drawn from the 
stellar evolution models of MazzitcUi & D'Antona 
(1986). In addition we have obtained (C. A. Tout, 
personal communication) the internal composi- 
tional profiles corresponding to the stellar evolu- 
tion models of Hurley et al (2000). We will use 
these two models as an indication of the variation 
to be expected due to uncertainties in the calcu- 
lation of interior compositions. 

How do the input models differ? Figure 24 
shows the mass fraction in carbon as a function 
of enclosed mass for a 0.6Mq white dwarf for the 
case of the Mazzitelli & D'Antona profile (solid 
line) and Hurley ct al profile (dashed line). We see 
that the default model has a higher carbon frac- 
tion in the inner parts and thus a correspondingly 
higher heat capacity (more ions per unit mass). 
Thus, we expect the Hurley models to cool a lit- 
tle bit faster. The cooling curves for two O.55M0 
models arc shown in Figure 25. The solid line is 
our default model and the dashed line is the model 
with the Hurley et al compositional profile. 

Figure 26 shows the fit for our default set of 
parameters, but now incorporating cooling mod- 
els using the Hurley et al internal profile. The 
best fit model is now {x,T) =(-0.95,12.4 Gyr) and 
the minimum Xdof ~ 0.8 is even lower than for 
the default model. The lower age limit is now 
10.4 Gyr and, with the lower x^, this limit will 
hold up after the marginalisation. The upper age 
limit moves to beyond > 15 Gyr (the 2a contour 
extends upwards beyond the left edge of the plot). 
Figure 27 shows comparison of the new best fit 
model with the observations. 

5.4. Initial-Final Mass Relation 

The extreme range of the variation of white 
dwarf mass in Figure 11 is given by our default 
model (which yields the lowest white dwarf mass 
at a given main sequence mass) and the pop II 
relation from Hurley et al (2000), which yields 
the largest white dwarf mass. If we replace equa- 
tion (8) with the latter relation but retain all 
other features of the default model, then we ob- 
tain the left-hand panel of Figure 28. Of course, 
a more self-consistent application is to apply the 
core composition from Hurley et al with the initial- 
final mass relation from the same models. The 



result is shown in right hand panel of Figure 28. 
There is little significant change in the overall fits, 
in either case. 

This is perhaps not too surprising for it is only 
for relatively young ages that the massive white 
dwarfs lie above our detection threshold and these 
ages are already excluded from our confidence re- 
gion. Figure 29 shows the results of our Monte- 
Carlo models for two different cases. On the left 
we show the white dwarf mass versus V magnitude 
for our best fit model, using equation (8). On the 
right we show the same thing, but for a model of 
age only 10 Gyr. A much larger fraction of the 
white dwarf population lies above our detcctabil- 
ity threshold and so we probe a much larger range 
of white dwarf masses. 

5.5. Other Group's Models 

As an extension of section § 5.3.2, we would 
like to compare with the models of other groups 
in the literature. We choose two sets of models 
from the recent literature - those of Chabrier et 
al (2000) and Salaris et al (2000). These have 
been chosen specifically because they incorporate 
not only up-to-date internal physics but also make 
an attempt to treat the outer boundary condi- 
tion properly using stellar atmosphere models. We 
know of two other calculations which meet this 
standard (Fontaine et al 2001; Prada Moroni & 
Straniero 2003) but which did not publish tables 
that could be used for this purpose. Figure 30 
shows the 0.6 Mq cooling curves for the three 
different groups. The most significant variation 
occurs at late times (T> 12 Gyr) where the dif- 
ferences are of order ±1 Gyr with respect to our 
models. 

In performing this comparison we encounter the 
problem that the published models are on a rather 
coarse grid in mass and thus not sufficient to prop- 
erly compare with the full analysis we publish 
above. To perform a fair comparison between dif- 
ferent models, we will first perform degraded fits 
to our models. We will assume all white dwarfs 
have mass 0.6 Mq (since all groups have published 
models for this canonical white dwarf mass), but 
we will keep the Ay = 1.38. Thus, strictly speak- 
ing the models will not fit the upper part of the 
main sequence correctly. However, most of the age 
discrimination in the models occurs at /> 26, so 
this is equivalent to assuming some modest evo- 
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lution in the initial-final mass relation such that 
M ~ O.6M0 at the magnitudes of interest. Fig- 
ure 31 shows the fit using our default cooling mod- 
els, but without enforcing the conditions of § 2.4. 
Comparing this fit with those below affords a di- 
rect comparison of the influence of other cooling 
models. The character of the solution is the same, 
with the lower age limit moving down to 10.4 Gyr. 
The best fit model has Xdof ~ 1-1' the model is 
still quite acceptable. 

The left panel of Figure 32 shows the fit ob- 
tained when we use the 0.6Mq cooling models of 
Chabrier et al (2000). Despite the inconsistency 
between the mass and the extinction, the Chabrier 
et al models also provide an acceptable fit, al- 
though with a marginal Xdof The allowed 
age range is somewhat lower (9.8-12.1 Gyr), not 
surprising given that these models are the fastest 
cooling of those compared in Figure 30. The right 
panel of Figure 32 shows the results obtained us- 
ing the models of Salaris et al. The best fit mod- 
els show an age range 14.4-16 Gyr, as appropriate 
to the slowest cooling models of those discussed. 
However, the minimum Xdof 2.1, indicating that 
these models do not provide a good fit. The fits are 
similarly poor with both the 0.55Mq and 0.6Mq 
models. 

It is encouraging that the fits obtained with 
our models and those of Chabrier et al have sim- 
ilar character and similar lower limits and that 
although the solutions obtained with the Salaris 
models arc markedly different they are also a very 
poor fit - demonstrating that the observations are 
now good enough to rule out some models. Per- 
haps in the end it is not surprising that the slow- 
est cooling models arc the most susceptible to con- 
straint - a greater fraction of the entire luminosity 
function lies above our observational threshold. 

5.6. Helium Atmospheres 

We have thus far considered only models with 
hydrogen atmospheres. However, observations of 
the Galactic disk show that a significant fraction 
of cool white dwarfs possess helium-dominated at- 
mospheres (Greenstein 1986; Bergeron et al 1990; 
Bergeron, Ruiz & Leggett 1997). The physics that 
determines this ratio is still somewhat uncertain 
and it remains an issue as to whether the appear- 
ance of atmospheric helium is a temporary phase 
in the life of a white dwarf - one in which the 



dwarf shows a helium atmosphere for a time, but 
whose cooling history is still similar to that of a 
dwarf with a hydrogen atmosphere boundary con- 
dition, or whether the cooling is truly that of a 
helium atmosphere white dwarf. This is an im- 
portant distinction because the low opacity and 
high photospheric pressures of a cool helium at- 
mosphere lead to much more rapid cooling. We 
will examine both cases, although most workers in 
the field favour the former. 

In the Galactic disk, the fraction of cool white 
dwarfs with helium atmospheres lies somewhere 
within the range 0.33-0.5. It is possible that ac- 
cretion from the ISM may be partially responsible 
for determining the balance between hydrogen and 
helium atmospheres (for a more detailed discus- 
sion see appendix A.) If so, the accretion history 
for white dwarfs in clusters may be significantly 
difii'erent and so we wish to examine a variety of 
possibilities. 

Thus we will treat the helium fraction as a free 
parameter. The proper procedure is then to evalu- 
ate the at each (x,T) value for a range of helium 
atmosphere fractions and then marginalise over 
the free parameter - adopting the best fit helium 
atmosphere fraction at each value of (x,T). The re- 
sulting fit is shown in Figure 33, scanning a range 
of hydrogen fractions from 0.0 to 1.0. The left- 
hand panel shows the resulting Hess fit contours. 
The solution looks very similar to the hydrogen 
case, shifting to only slightly lower ages. The best 
fit model is now at (a;, r)=(-1.2,12.4 Gyr), and 
is found for a model in which 60% of the white 
dwarfs cool with hydrogen atmosphers and 40% 
with helium atmospheres. This is consistent with 
the results found in paper I. in which there was lit- 
tle change in the model as long as more than 50% 
of the white dwarfs have hydrogen atmospheres. 
This is primarily because the most striking feature 
of the observed luminosity function is the sharp 
rise in the numbers at / > 26. As we have de- 
scribed in § 4.4, this jump is a consequence of the 
slowing in cooling due to the change in the outer 
boundary condition that results when the photo- 
spheric opacity drops due to the recombination of 
atmospheric hydrogen. When the atmosphere con- 
sists of helium instead, the evolution of the outer 
boundary condition is different. In particular, re- 
combination occurs at higher efisctive tempera- 
tures. Consequently, we require a model with a 



17 



significant fraction of liydrogcn atmospheres in or- 
der to matcli tiie observations. 

If some white dwarfs change their atmospheric 
appearance as they age, then the cooling ages may 
still be more representative of a hydrogen atmo- 
sphere boundary condition even if they do appear 
to show a helium atmosphere at some point in 
their lives. In this case, a more accurate reflec- 
tion of reality is to simply assume all white dwarfs 
cool according to a hydrogen atmosphere bound- 
ary condition, but to assume that some fraction 
show colours more characteristic of helium atmo- 
spheres. Figure 34 shows the results of a parame- 
ter scan if we marginalise over the fraction of stars 
with helium atmosphere colours only (i.e. we as- 
sume they all still cool with hydrogen atmosphere 
boundary conditions). Again, the results do not 
differ significantly from the default model. 

5.7. Sensitivity to atmosphere models 

We have also calculated fits using the atmo- 
sphere models of Hansen (1998, 1999). Figure 35 
shows the fit for our default model. The best 
fit model moves slightly downwards to {x,T) =(- 
1.2,12.0 Gyr) and the confidence limits appear to 
tighten somewhat. However the minimum Xdo/ ^ 
2.1, so that this is not as good a fit as with the 
Bergeron atmospheres. 

5.8. Binary Fraction 

We inferred in § 4 that our white dwarf sam- 
ple contains a binary fraction of ~ 10% or more. 
A non-zero binary fraction may also have an in- 
fluence on our results because it makes stars ap- 
pear anomalously bright at fixed colour. We ad- 
dress this by allowing some given fraction of white 
dwarfs in our sample to be unresolved binaries. 
The binaries are assumed to be randomly sampled 
from the isolated white dwarf distribution i.e. we 
do not take into account any influence of binarity 
on the overall stellar evolution of the white dwarf. 
Note, we only consider double degenerate bina- 
ries - even the lowest main sequence star is bright 
enough to dominate a white dwarf companion suf- 
ficiently to move the system beyond our colour 
cuts. The solution for a 10% binary fraction does 
not change significantly, with a best fit {x,T) =(- 
1,12.4 Gyr) and a lower limit of 10.8 Gyr. 



5.9. Proper Motion Cutoff 

Figure 3 shows that our proper motion accu- 
racy gets worse as the stars get fainter. Eventu- 
ally, the scatter of the white dwarfs becomes com- 
parable to the cut wc assign to separate cluster 
members from background and thus we must ex- 
amine how sensitive our results are to the proper 
motion cut. All the preceding calculations have 
used the proper motion cut A in Figure 3. In 
Figure 36 we show the results obtained if we use 
the more stringent proper motion cut designated B 
for our default model. Note that we do the com- 
parison self-consistently in that we recompute all 
the artificial star tests with the new proper motion 
limit. The age range does not change significantly, 
although the solution is shifted to slightly larger 
values of x. Some expansion of the confidence con- 
tours is expected as we have reduced the number 
of observed stars by 16%. 

5.10. Distance Uncertainties 

The uncertainty in the distance to globular 
clusters is the principal source of error for the main 
sequence fitting method. As outlined in § 2.4, dis- 
tance errors and extinction are intimately tied in 
our analysis method. Keeping the white dwarf 
mass at the top of the cooling sequence fixed, a 
change in the absolute distance modulus is re- 
flected in the extinction. A 5% reduction in the 
distance for our default model reduces the extinc- 
tion to Ay = 1.12. Recall that the counterin- 
tuitive fact that extinction changes in the same 
sense as the absolute distance modulus is because 
of the change in the reddening dominates over the 
change in the extinction in matching the observed 
sequence (§ 2.4). 

Thus, to investigate the effect of distance un- 
certainty we allow for a range in distance and ex- 
tinction, but require both to meet observational 
bounds. We consider the true distance modulus to 
vary in the range /z = 11. 18 ±0.18 (so the distance 
uncertainty is less for M4 than for most other clus- 
ters, because its proximity to the sun allows astro- 
metric measures to be made) and consider a loose 
bound on the extinction of 1.1 < Ay < 1.5. We 
have allowed errors twice as large as inferred from 
the subdwarf match, to account for possible dif- 
ferences between white dwarf and main sequence 
spectral shapes. We investigate a series of masses 
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0.5 < M < 0.6Mq and for each mass examine 
the range of /x and Ay that satisfies these con- 
straints. This is shown in Figure 37 - we see that 
the extinction limits rule out the lower part of the 
distance range for O.SM© and the upper part of 
the distance range for Q.QMq. Masses larger than 
O.65M0 arc ruled out by this procedure. 

Figure 38 shows the result of marginalising our 
model fits over the allowed distance and extinc- 
tion range for the full range of upper cooling se- 
quence white dwarf masses Q.^Mq-Q.^Mq. Also 
shown is the 2(j contour for the allowed range of 
distance and extinction but restricting the mass 
to O.55M0. The lower age limit thus obtained is 
11.1 Gyr. The fact that this is larger than that 
obtained for just the O.SSMq model is a reflection 
of the fact that the minimum is found with 
M=O.6M0 and = H-l- The location of the 
minimum is again similar to that of prior cases 
(a;,T)=(-l.l,12.1Gyr). 

6. Discussion 

We have presented a rather bewildering ar- 
ray of possible systematic uncertainties which can 
influence the model of the cluster white dwarf 
sequence. To properly assess the global sys- 
tematic error on our age determination, we now 
marginalise over all the model variations. We in- 
clude all the variations in physical inputs on our 
model, including core compositions, atmospheric 
models and binary fraction. We do not include 
the fits obtained using other models from the lit- 
erature, as we wish to express the quantifiable 
model systematics within our model alone. The 
result is shown in Figure 39. The global mini- 
mum is found at (x,T) = (-1.2,12.1 Gyr), with a 2a 
lower limit of 10.3 Gyr. A summary of the various 
model fits is given in Table 1. Figure 40 shows the 
Hess fit for our best fit model. Although we have 
not used the luminosity function fitting method 
to locate this minimum, the resulting luminosity 
functions, shown in Figure 41 arc also good fits, 
with = 7.9 for V and = 3.8 for 1, where each 
has 5 (7-2) degrees of freedom. 

The overall character of the solution remains 
unchanged throughout the variation of the sys- 
tematic uncertainties. A strong covariance be- 
tween age and the parameter x is evident, in the 
sense that larger ages require smaller x, with ages 



~ 14 Gyr requiring x ^ —2. A shallower 'slope' 
of X ~ implies a best fit age ~ 11 Gyr The 
largest x solutions have x ~ 0.5, with an age 
T ^ 11.3 Gyr. A comparison between the white 
dwarf and main sequence number counts yields 
a similar x-T covariance (§ 5.2), because larger 
ages imply a greater fraction of the total white 
dwarf population has cooled beyond our detection 
threshold, requiring a shallower progenitor mass 
function to match the observed counts. We should 
note that the full calculation does not use any in- 
formation about the main sequence number counts 
{x is used only to choose the probability distri- 
bution of progenitor masses in the Monte Carlo 
model), so that the fact that we obtain the same 
trend from two different lines of calculation is a 
sign of consistency. We repeat again our caution- 
ary note from § 5.2, that x should be considered 
only as a parameter in the probability distribution 
of white dwarf progenitors. The relationship be- 
tween X and the true mass function slope for the 
global cluster population is affected by the dynam- 
ical processes in the cluster. As we have shown, 
Michie-King models for this cluster indicate that 
the true mass function has a slope x > 1. 

Our result is thus a best fit age of 12.1 Gyr, 
with a 95% lower bound of 10.3 Gyr for the age of 
M4. This is remarkably similar (12.6 Gyr best fit, 
with 95% lower bound 10.4 Gyr) to that quoted by 
Krauss & Chaboyer (2003) for the cluster system 
as a whole. There is some evidence that there is 
some dispersion in age amongst the clusters, par- 
ticularly for the somewhat more metal-rich group 
(Stetson, VandenBerg & Bolte 1996; Rosenberg et 
al 1999) which includes M4. However, our meth- 
ods are not yet accurate enough to place a strong 
constraint on this question. It is also encourag- 
ing that solutions of similar character are obtained 
with the completely independent models published 
by Chabricr et al (2000). Even the fact that wc 
are unable to achieve a good fit with the models 
of Salaris et al (2000) is encouraging - it suggests 
that our data and methods arc now robust enough 
to begin discriminating between models. This is a 
necessary step in the direction of making the white 
dwarf cooling method an alternative age discrimi- 
nator. 

Our age determination for M4 using the white 
dwarf cooling sequence allows us to perform a 
direct comparison between the age of the solar 
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neighbourhood and the globular cluster system, 
using the same method (paper I). Our best fit 
model has not changed significantly since paper I, 
but the larger error bar due to model uncertain- 
ties has brought the lower limit down to 10 Gyr. 
This is still larger than the 7.3 ±1.5 Gyr {2a 
range) we infer for the solar neighbourhood. Al- 
though the best fit models for the two popula- 
tions support the considerable delay predicted by 
some models for thin disk formation (e.g. Burk- 
ert, Hensler & Truran 1992), our 2a limits are now 
consistent with a delay ~ 1 Gyr. It is also worth 
noting the dominant systematics are different in 
the two cases. In the case of the Galactic disk, 
the faintest white dwarfs are dominated by he- 
lium atmospheres, which appear to be a minority 
constituent in the globular cluster model (§ 5.6). 
Figure 42 shows the sample of disk white dwarfs 
from Liebert, Dahn & Monet (1988), when we ap- 
ply our preferred value of extinction and distance 
to place them on the M4 cooling sequence. The 
turnover in the disk luminosity function occurs at 
Ml ~ 14.3 (which corresponds to / ~ 26.4 for our 
default parameters), well above the faintest white 
dwarfs in our sample a clear indication of the 
age difference between the cluster and the Galac- 
tic disk. 

Our description of Galactic chronology is con- 
sistent with other estimates as well. The ages of 
solar and super-solar metallicity field stars with 
Hipparcos parallaxes appear to be ~ 7-8 Gyr (Liu 
& Chaboyer 2000; Sandage, Lubin & VandenBerg 
2003) although some metal poor stars with larger 
ages (^ 10 Gyr) may be kinematically part of the 
thin disk. 

To conclude, in this paper we have performed 
a detailed analysis fitting models for the cooling 
history of the cluster white dwarfs to our obser- 
vations. We emphasize again that we fit models 
to the full luminosity function - we do not rely on 
localised features such as a putative truncation in 
the luminosity function. We have demonstrated 
that not only is it possible to get excellent fits to 
the observations, but it is also possible to rule out 
significant portions of parameter space. This lat- 
ter point is important as it demonstrates that the 
method has some power to discriminate between 
models. We stress this because in the past the 
comparison with white dwarf models and observa- 
tions has been rather more qualitative than quan- 



titative. This was in part due to the poor statistics 
at the faint end of the disk luminosity function and 
in part because the uncertainties injected into the 
model by an extended star formation history. The 
quality of observations has now reached the level 
where more detailed and statistically rigorous the- 
oretical analyses are possible. This is a necessary 
step in making the white dwarf cooling sequence 
a useful tool for age determination. 
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A. Chemical Evolution of White Dwarfs 



It is an established observational fact that a significant fraction of cool white dwarfs have helinm-rich 
atmospheres. Indeed, there is much evidence to suggest that the chemical composition of an individual 
white dwarf evolves as it cools (Greenstein 1986; Bergeron et al 1990; Bergeron, Ruiz & Leggett 1997). The 
first clement is that the fraction of helium atmosphere dwarfs increases at cooler temperatures. Thus, some 
hydrogen atmosphere dwarfs become helium atmosphere dwarfs. This is thought to be due to the growth of 
the surface convection zone, which can dredge up helium from deeper layers and mix it with a thin surface 
hydrogen layer. However, it appears that at lower temperatures the composition reverts back to a hydrogen 
atmosphere, at least temporarily. The reasons for this evolution are still not well understood, but there is 
general agreement that it has to do with the combination of the depth of the convection zone and the amount 
of hydrogen on the surface. Indeed, the accretion of hydrogen from the interstellar medium may play an 
important role. 

If accretion from the ISM is important, then we might reasonably expect differences in the chemical 
evolution of population I and population II white dwarfs, as the former spend more time in the high density 
ISM of the Galactic disk and move will lower velocities. On the other hand, white dwarfs in clusters may 
have even higher accretion rates if the intracluster medium is enhanced as a result of cluster mass loss. 

We may estimate the effect of the cluster environment by converting our present-day global white dwarf 
birthrate (~ 10~^yr~^ - see § 4.1) into a mass deposition rate into the cluster ISM (~ 3 x 10~'^ M(7,.yr~^). 
If this mass builds up in the cluster potential in between disk passages (which most likely result in periodic 
stripping of the cluster gas) then we expect roughly ~ 30Mq to have built up in the cluster ISM, leading to 
a mean density n ^ 3cm~^ (if we use the half-mass radius ~ 3 pc). This is actually higher than the solar 
neighbourhood density. Furthermore, the stellar velocity dispersion is lower in the cluster as well, so that 
the accretion rate is 



approximately 600 times larger than for a corresponding white dwarf in the solar neighbourhood. Note 
that the resulting accretion Imninosity in this case is 4 x 1O~^L0 (an issue first raised by Lin, personal 
communication), which translates to a magnitude V ^ 28.8 - i.e. right near the limit of our detection. 

However, this is almost certainly an over-estimate. The fact that the clusters do not exhibit a substantial 
intracluster medium has been known for many years (Hesser & Shawl 1977; Gnapp et al 1996; Hopwood et al 
1999). Only recently intracluster gas actually been conclusively detected (Freire et al 2001). In this case, the 
cluster is 47 Tuc, and the density ~ 0.07cm^^ corresponds to only ~ O.IMq in the inner 2.5 pc of the cluster. 
Although this is determined from pulsar dispersion measures and hence only a measure of the ionized gas, 
the UV flux resulting from the high stellar density is expected to ionize most of the cluster gas, so we will 
adopt a more conservative measure of the M4 intracluster density as O.lcm"^. Furthermore, if the gas is 
predominantly ionized, the thermal velocity (~ 9km. s~^) is larger than the stellar velocity dispersion, which 
further reduces the expected accretion rate, to 2 x 10~^^MQ.yr^^ (and reduces the accretion luminosity 
to well below our detection limit). This rate is still a factor ~ 3 larger than expected for the Galactic disk. 

Thus we conclude that, if accretion has a significant influence on the atmospheric chemical evolution of 
white dwarfs, then cluster white dwarfs probably accrete more hydrogen than even disk white dwarfs. Thus, 
it is quite possible that the temperature range over which disk dwarfs show a dearth of hydrogen atmospheres 
(Bergeron, Ruiz & Leggett 1997) may be smaller for the cluster white dwarfs and that the overall fraction 
of dwarfs with hydrogen atmospheres may be larger for the cluster. 

This should also introduce a note of caution if we wish to use the cluster white dwarfs to make any 
assertions about putative populations of white dwarfs in the Galactic halo. While the two populations 
do indeed share similar progenitor histories, their subsequent accretion histories are likely to be markedly 
diflferent. Cluster white dwarfs most likely enjoy an abundance of accreted hydrogen even in the face of the 
constant removal of gas from the cluster, whereas true halo white dwarfs spend most of their lives in the 
very underdense regions of the Galactic halo. 




(Al) 
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B. The Gaussian Nature of the Errors 

An implicit assumption in the use of the statistic to determine confidence intervals is that the distri- 
bution of errors is Gaussian^^(e.g. Press et al 1992). Since we are dealing with number counts of stars one 
might naturally assume this is true (or at least that the errors are Poissonian, if not Gaussian). However, 
we note that much of the spread in the white dwarf cooling sequence is the result of photometric errors. We 
wish to examine whether the Poissonian nature of the number counts is still to be expected when we include 
the artificial star tests in our calculation. 

The variable whose probability distribution we need to examine is Ni, the number of stars counted in a 
given bin, i. Thus, we sample our model distribution a fixed number of times (A^o) and then run it through 
the probability distributions that result from our artificial star tests. We then examine the distribution of N^. 
The true value A^o is unknown, so we fix it such that it yields a prescribed mean value of Ni . Figure 43 shows 
the resulting cumulative distribution (solid histograms) of Ni for two of our faintest bins, 26.5 < I < 27.0, 
1.45 <V - I < 1.75 and 26.5 < I < 27.0, 1.75 <V - I < 2.05, when the input white dwarf was located at 
the center of the first bin. The mean values are arbitrary, but of the same order as the true number counts 
and specifically chosen so that the two distributions do not coincide too closely. The dotted lines indicate 
the cumulative Poisson distributions for the respective mean values of each bin. We see that the errors are, 
if anything, slightly better than Poissonian! 

C. An 'alternative analysis' 

Recently, de Marchi et al (2003) have published a re-analysis of our data which fails to reproduce our 
result. We have severe reservations about this work on both observational and theoretical grounds. 

On the observational side, de Marchi et al find that they cannot achieve field-cluster separation at the 
depths we do. In fact, their limit is almost a magnitude brighter than ours. The detailed reasons for this 
difference are discussed in Richer et al (2004b). In short, de Marchi et al use image-association stacks 
generated by the Canadian Astronomy Data Center and the ESO Space Telescope Coordinating Facility 
(Micol & Durand 2002). The resulting images of faint sources are significantly more diffuse than those 
determined by our method, and so the positional accuracy of the de Marchi et al images is cleary degraded. 
However, the principal fault with the analysis of de Marchi et al is their refusal to use the recentering 
algorithm in ALLSTAR, which severely limits their ability to perform proper motion separations. They 
claim that this is prone to measuring noise spikes as cluster members, but it is easy to demonstrate by a 
simple experiment that this is an unfounded concern. By using the original finding list but randomising the 
stellar positions. Richer et al (2004b) show that very few of these spurious 'stars' pass the photometric cuts 
we apply and those that do are rejected by the cluster proper motion cut. 

Given the fiawed data analysis of de Marchi et al there is little point in a detailed critique of their 
theoretical analysis. However, a few points are worth noting. As we discuss in § 2.4, it is important for 
internal consistency to ensure that the model white dwarf sequence actually overlaps the observations for a 
given choice of mass, distance and extinction, de Marchi et al quote a single choice for extinction (Harris 
1996) which requires a white dwarf mass ~ O.6M0 at the bright end, yet they quote a range of white dwarf 
masses. So, while some of their models are clearly consistent at the bright end, it isn't clear that they all are. 
It is also worth noting that some of the alternative extinction and distances they mention (e.g. Djorgovski 
et al 1993) can be ruled out as they would require unphysically low mass (< O.5M0) white dwarfs. 

The statement that their photometric errors were drawn from a Gaussian raises concern. As can be seen 
from Figure 6 the dispersion at fixed magnitude can be different for stars recovered brighter and fainter than 
the input magnitude. Furthermore, the wings of the distribution are often larger than Gaussian - which 
is important to include in the analysis. We note here we are talking about the distribution of photometric 

'The determination of the best fit model (the minimum of x^) is not affected by this, but it is necessary to assign statistical 
significance to a value of Ax^. 
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errors, not the distribution of number counts in a given bin, which we have demonstrated (§ B) is indeed 
well approximated by a Gaussian. 

We have demonstrated in this paper that fitting to the full Hess diagram is a much better method for 
model fitting than comparing to one-dimensional luminosity functions, de Marchi ct al state a preference for 
the I LF that may bias their results as our experience suggests fitting to the I LF gives consistently lower 
ages than either the V LF or the full Hess fit. It is interesting to note that they furthermore claim that they 
require an 'unphysical' value of a = (corresponding to a; = —1) in order to match our luminosity function. 
As we have discussed in § 5.2, the assumption that a or a: is indeed the true mass function slope is a dangerous 
one, as this cluster clearly shows significant mass segregation. As a result, the white dwarf population in 
our field is a product not only of stellar evolution but cluster dynamical evolution as well. The quantities a 
or X should be viewed as parameters in a probability distribution foremost, and their interpretation as true 
physical values should be treated with great caution. 

We close with a brief philosophical point. De Marchi et al assert that they would consider the clear location 

of the luminosity function peak as a mandatory requirement for determining the cluster age. If one is going 
to determine a cluster age, it will be the result of a model fit and we fail to see why one would believe an age 
determined from fitting a model to the cutoff if one does not believe that the model also accurately describes 
the luminosity function above the cutoff! While the determination of a Galactic disk age has generally been 
performed by this method, that was largely because a proper and rigorous statistical treatment of the data 
is hampered in that case by poor statistics and the convolution of white dwarf cooling with an uncertain 
star formation rate. We have demonstrated in this paper that the current level of observations in M4 is 
sufficient to allow us for the first time to perform proper, statistically rigorous, model comparisons. The 
demonstration of this is that we can not only find models that fit the data but furthermore can rule out 
significant portions of parameter space. We believe it is time to move beyond the hidebound attitude that 
only the faintest white dwarf in a given population can be used for age discrimination. 
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Fig. 1. — The proper motion displacements dx and dy (in units of HST pixels) are shown for the 6 year 
baseline spanned by our observations. The co-ordinate system is centered on the bright cluster stars. All 
objects in the field are shown on this diagram. A simple proper motion cut of 0.5 pixels separates cluster 
members (filled circles) from background (open circles). 

Fig. 2. — The colour magnitude diagram for all detected sources is shown. Three principal structures are 
evident - the cluster main sequence at the upper right, the inner halo (i.e. background) main sequence in 
the middle, and the cluster white dwarf sequence at the lower left. The dashed line indicates a generous 
cut to define the region of the CMD that will contain cluster white dwarf candidates. The proper motion 
displacements for all sources below this line are shown in figure 3. 

Fig. 3. Here we show the proper motion displacements of all the 'cluster white dwarf candidates' of 
Figure 2. The horizontal axis /x is the total proper motion displacement relative to the frame defined by 
the bright cluster main sequence. It clearly splits into two groups, cluster white dwarfs (on the left) and 
background objects. Even at faint magnitudes, the concentration of displacements to the left indicates that 
most objects are cluster white dwarfs. To investigate the infiuence on our results of which proper motion 
cut we use, we define two different proper motion selection criteria, marked as A and B. Another point to 
note is that considerations of completeness will restrict our attention to white dwarfs with I < 27. 

Fig. 4. — The resulting proper motion-selected white dwarf cooling sequence for Messier 4 is shown. Also 
shown is the reddening vector - thus a change in the extinction moves stars along the sequence much more 
than across it. 

Fig. 5. — The three different curves indicate the fraction of stars recovered as a function of magnitude 
for all three WFPC2 chips. This fraction includes stars recovered at any magnitude but does require that 
the position be within 0.5 pixels of the true position i.e. lie within the proper motion cutoff. The 50% 

completeness limit is located at ^ 29 for all three chips. 

Fig. 6. — The points represent the la bounds on the distribution of recovered magnitudes (relative to input 
magnitude) as a function of input magnitude, for the V band (upper panel) and I band (lower panel). The 
solid points represent the tail to brighter magnitudes (so 16% of artificial stars were recovered at a magnitude 
brighter than —a above the input magnitude) and the open circles the tail to fainter magnitudes (so 16% 
lie fainter than mag-|-cr). The dotted line indicates the 50% completeness limit in both bands. Recall again 
that all artificial stars counted also have to pass the proper motion cut. 

Fig. 7. — The upper part of the cooling sequence (solid points) is shown compared to the theoretical model 
fit (dashed line) for the case of a O.55M0 model. The line can be adjusted not only vertically (by distance 
and change in mass) but also with a horizontal component (due to the changes in reddening associated with 
changes in extinction). 

Fig. 8. — The solid curve shown is the cooling sequence for a white dwarf of fixed radius. The value chosen 
in this case was appropriate to a cool O.55M0 dwarf and thus the extinction is inferred to be Ay = 1.39 to 
fit the upper part of the cooling sequence. 

Fig. 9. — The three models shown are identical O.55M0 cooling models except for the different hydrogen and 
helium layer masses. Our default values are for mass fractions qh = 10~^ and qne = 10~^. These are the 
most commonly used values and also at the 'thick' end of the expected spectrum. More massive hydrogen 
layers are unlikely (it would result in significant nuclear burning i.e. the star would become a red giant 
again) but they could potentially be thinner. We show the cooling curve for a qh = 10~^ as well (dotted 
curve). The effects of a thinner helium layer are also shown (dashed curve). 

Fig. 10. — The solid line represents the duration of the pre-white dwarf evolutionary phases as a function of 
Zero-age main sequence mass for stars of low metallicity appropriate to M4. The dotted line indicates the 
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corresponding value for a solar metallcity sequence. Both sets of models are taken from the results of Hurley 

et al (2000). 

Fig. 11. — The solid line is model A from Wood (1992), which is the model for default relation (8). The 
dotted (Zq) and dashed (O.IZ0) are relations from Hurley et al (2000) which indicate the range of variation 
expected from changes in metallicity between population I and population 11 white dwarfs. The solid points 
are the Z = O.OSZ© models of Dominguez et al (2000), while the asterisks represent the empirical relation 
from Weidemann (2000). 

Fig. 12. — The left-hand panel shows the cooling curve for a 0.55Mq model. The right hand panel shows 
the same cooling sequence in the colour-magnitude diagram, compared to the observed cooling sequence. 
The arrows indicate the corresponding positions of two important epochs. The first is the onset of core 
crystallisation, which contributes latent heat. The second is the point at which the effective temperature 
reaches 5000 K. This corresponds to a drop in photospheric hydrogen opacity and a change in the boundary 
condition of the cooling models. The result is a flattening of the cooling curve which causes a corresponding 
jump in the observed hmiinosity function. 

Fig. 13. — The different curves indicate the cooling curves for white dwarfs of different mass. At early times 
more massive white dwarfs are brighter because their larger mass corresponds to a larger heat capacity. 
However, the higher central densities means that more massive white dwarfs crystallize sooner and at late 
times cool faster because they enter the Debye cooling regime earlier. 

Fig. 14. — The left-hand panel shows the upper cooling sequence for the annulus fro 1-2 core radii of the 
cluster M4. The middle and right-hand panels show the same for the ranges 2-3 and 3-5 core radii. The 
horizontal dotted lines bracket the region in which we determine the binary fraction. The dashed line is 
the designated upper envelope of single white dwarfs. Our binary fraction is determined by counting what 
fraction of the white dwarfs in this magnitude range lie to the right of the dashed line. The error bars in the 
left and middle panels indicate the photometric uncertainties at given places along the cooling sequence. 

Fig. 15. — The boxes delineate the grid elements for the fits to the full colour distribution (Hess fits). The 

upper cooling sequence is defined as a single grid point since there is little age information in the detailed 
distribution at these bright magnitudes. The dashed lines indicate the 50% completeness limits for V and I. 

Fig. 16. — The left-hand panel shows the 68% (heavy solid contour) and 95% (thin solid contour) confidence 
intervals that result from the Hess diagram fits to our default model. The right hand panel shows the 
corresponding results if we fit only to the V band luminosity function (solid contour) or the I band luminosity 
function (dotted contour). 

Fig. 17. — The upper left-hand panel shows the number counts in the large, bright white dwarf bin. The 
other three bins compare the colour distribution in the labelled magnitude ranges with the observations. 
The models match both the ratio of bright and faint white dwarfs as well as the colour distribution of the 
faint white dwarfs. 

Fig. 18. — The upper panel shows the V-band luminosity function and the bottom panel shows the corre- 
sponding I LF. The dashed lines in each case enclose the region where the fit was performed. Fainter bins 
were ignored because of poor completeness, while bins brighter than this region were summed together and 
included as a single bin in the fit. 

Fig. 19. — This fit is for a young (10 Gyr) solution and clearly the model contains too many faint white 
dwarfs relative to the bright end. Furthermore, the colour distribution is too skewed towards the blue at the 
faint end. 

Fig. 20. — In this case we have the opposite problem to that in Figure 19 - namely that there are too many 
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bright white dwarfs relative to the faint white dwarfs. The colour distribution is also flatter than that shown 

by the observations. 

Fig. 21. — The solid contours show the 1 and 2 a bounds from the Hess fits for our default model. The 
dashed lines indicate the relation between x and T we obtain if we match the white dwarf and main sequence 
number counts, as described in § 5.2, for three different assumed break masses. 

Fig. 22. — The model for a thin hydrogen envelope gives a similar best fit age but the lower age limit 
decreases. The solid cm-ves show the Icr and 2<t contours for this fit, compared with the 2cr contour from 

Figure 16 (dotted line). 

Fig. 23. — A thinner helium layer increases the age bounds slightly. Once again, the solid contours indicate 
the 1 and 2 a intervals for the model under discussion and the dotted contour shows the 2a interval for our 
default model, included as a point of reference. 

Fig. 24. — The solid line shows the profile of carbon fraction Xq as a function of mass enclosed for our 
default O.6M0 model. The dashed line shows the profile we obtain from the Z=0.001 Hurley et al models. 

Fig. 25. As in Figure 24, the solid line indicates our standard model, while the dashed line is for the 
model using the Hurley et al (2000) inputs. The greater oxygen content of the latter model results in earlier 
crystallisation and thus faster cooling at late times. 

Fig. 26. — The Hurley et al core composition shifts the 2(t confidence region (thin solid contour) to slightly 
lower values than the equivalent region for the default models (dotted contour). 

Fig. 27. — Here we compare the model for the case (a;, r)=(-0.95,12.1 Gyr) with the observations. The 
colour distributions and the relative normalisations of the different bins all provide excellent fits. 

Fig. 28. — The left-hand panel shows the results of combining the Hurley et al initial-final mass relation 
with our default model. The dotted line indicates the 2cr contours of the default model, demonstrating that 
there is very little difference. The right hand panel shows the same thing but now using the core composition 
from Hurley et al as well. 

Fig. 29. — The left hand panel shows the I magnitude versus white dwarf mass for a model of age 12.4 Gyr. 

The right hand panel shows the same for a model of 10 Gyr. In both cases, the main sequence-white dwarf 
mass relation was defined as in equation (8). Each of the two panels contains 300 model stars. The scatter 
in the diagrams is the result of photometric errors derived from our artificial star tests. The dotted line 
indicates the 50% completeness limit used in our Hess fits. These diagrams illustrate that, in the case of 
older solutions, the more massive white dwarfs have cooled beyond the detection limit and we probe only a 
limited range of white dwarf masses. Consequently our results are not very sensitive to the exact form of 
the main sequence-white dwarf relation. 

Fig. 30. The solid line is our default cooling model and the dashed line is the same but with the core 
composition of Hurley et al (2000). The the filled points are from Chabrier et al (2000) and the open points 
are from Salaris et al (2000). At late times, these other two solutions bracket ours, with the Chabrier models 
cooling the fastest and the Salaris models the slowest. The latter would cool even more slowly if we included 
their calculated delay due to phase separation. 

Fig. 31. This shows the solution obtained, using our models but under the approximation we have to 
adopt to calculate the contours shown in Figure 32. In particular the assumption of O.6M0 but using our 
extinction of Ay = 1.39 means that the upper part of the cooling sequence is not accurately matched. The 
contours are similar to before, but the best-fit is slightly worse. 

Fig. 32. — The left hand panel shows the fit for the Chabrier et al (2000) models. The fit is acceptable. 
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despite the necessary approximation of 0.6Mq with Ay = 1.39. Again the solution shows the famihar 
covariance between x and T. The 2a hmit in this case extends down to ages of 9.8 Gyr. The right hand 
panel shows the corresponding fit for the Salaris models. The for these fits does not correspond to a good 
model. 

Fig. 33. — The solid contours shows the usual confidence regions based on the fit to the observed Hess 
diagram, now marginalised over Helium atmosphere fraction at each point. The dotted contour shows 
the 2(7 contour from the default model. The small difference indicates that the uncertainty in the helium 
atmosphere fraction does not have a big effect on the determination of the best model. 

Fig. 34. — This figure is the same as Figure 33, except for the fact that now the helium models still cool with 

outer boundary conditions determined by hydrogen atmospheres, although their colours are calculated with 
helium atmospheres. This is designed to mimic the case when the helium atmosphere appearance is only a 
temporary phenomenon in the life of a white dwarf. Again, the influence of the atmospheric uncertainties is 
small. 

Fig. 35. — Using the atmosphere models of Hansen (1998) we find a similar solution to before, but with 
a worse x^- The higher minimum chi"^ is responsible for the narrower confidence intervals, since the latter 
correspond to fixed Ax^. 

Fig. 36. — The solid curves indicate the confidence intervals if we evaluate our default model but now using 
the proper motion cut labelled B in Figure 3. The dotted curve is once again the 2a limit for the default 
model using proper motion cut A. The age range is not significantly affected, but the range of x values does 

now extend to larger values. 

Fig. 37. — The dotted lines delineate the bounds on the distance and extinction to M4. The solid lines 
indicate the relation between n and Ay that is required to fit the upper part of the cooling sequence, for 
each stated white dwarf mass. The cross indicates the location of the values preferred by Richer et al (1997). 
This lies close to our current default model (/i = 11.18, Ay ~ 1.39). 

Fig. 38. — The solid contours indicate the usual 1 and 2a bounds obtained by marginalising over the range 
of masses 0.5Mq-0.6Mq and the range of distance and extinction outlined in Figure 37. The dotted contour 
is different from our usual reference model, and shows the 2a range if we restrict our attention to O.SSM© 
but still allow the distance/extinction to vary. 

Fig. 39. The solid lines indicate the la (heavy line) and 2a (thin line) confidence intervals when we 
marginalise over our various systematic model uncertainties. The two asterisks indicate the locations of the 
two almost equal minima in the x^ surface. 

Fig. 40. — Our best fit model clearly meets all the various requirements imposed by the observations. The 
overall ratio of bright white dwarfs to faint white dwarfs is matched while the colour distribution of the faint 
white dwarfs is also well described. 

Fig. 41. — The upper panel shows the V LF while the lower panel shows the I LF. For both measures the 
rise in the luminosity function is well matched as is the plateau before the number counts drop off due to 
the incompleteness. Thus, although we have used the Hess diagram to find our best model, we clearly get a 
good fit to both of the luminosity functions as well. 

Fig. 42. — The left hand panel shows the M4 white dwarf sequence, with our analysis grid superimposed. 
The right hand panel shows the same grid but with the solar neighbourhood white dwarfs (Liebert et al 
1988) plotted at the location of M4. We have used the photometry of Leggett, Ruiz & Bergeron (1998) 
and an extinction of Ay — 1.39. The filled circles indicate dwarfs with hydrogen atmospheres and the open 
circles indicate those with helium atmospheres. Recall that the faintest 8 white dwarfs fall beyond the peak 
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of the disk luminosity function. The fact that the brighter of the disk white dwarfs he at the edge of our 
upper bin is an indication of the mass difference between disk white dwarfs and the brightest cluster white 
dwarfs. 

Fig. 43. — The solid histograms indicate the probability distribution of number counts for a fixed number of 
input stars. The stars were input with 1=26.7, V-I=1.75. The two bins cover the colours as indicated and 
26.5 < I < 27.0. The dotted lines are the incomplete T functions Q(N,17) and Q(N,26), corresponding to 
the appropriate cumulative Poisson distributions. 
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Table 1 
M4 Parameter Scan 



Mq/Av 


Cooling 


Qh Iqhp 


Atmosphere 


Best Fit 




Age (95%) 












(Gyr) 


min 


Gyr 


O.55M0/I, 


.385 


Hansen (default) 


10-4/10-2 


Bergeron H 


12.4 


0.97 


10.9-14.0 


0.55Mq/1, 


.385 


Hansen (default) 


io-Vio-2 

10-4/10-3 

10-4/10-2 
10-4/10-2 


Bergeron H 


12.4 


1.14 


9.8-13.6 


O.55M0/I, 


.385 


Hansen (default) 


Bergeron H 


11.7 


1.13 


11.0-14.2 


O.55M0/I, 


.385 


Hansen (Hurley) 


Bergeron H 


12.1 


0.77 


> 10.5 


O.6OM0/I, 


.385 


Hansen (default) 


Bergeron H 


12.4 


1.09 


10.4-14.5 


O.6OM0/I, 


.385 


Chabrier 


10-4/10-2 


Bergeron H 


10.2 


1.7 


9.8-12.1 


O.6OM0/I, 


.385 


Salaris 


10-4/10-2 
10-4/10-2 


Bergeron H 


15.5 


2.1 


15.0-16.3 


O.55M0/I, 


.385 


S alar is 


Bergeron H 


15.5 


2.1 


14.5-16.0 


O.55M0/I, 


.385 


Hansen (default) 


fn 10-4/10-2 


Bergeron /jj H 


12.4 


0.98 


10.5-14.0 








l-fH 0/10-2 
10-4/10-2 


& 1 - ,fo He 








O.55M0/I, 


.385 


Hansen (default) 


Bergeron Jh H 
& 1 - /h He 


14.6 


0.9 


> 11.6 


O.55M0/I, 


.385 


Hansen (default) 
Hurley IFMR 


10-4/10-2 


Bergeron H 


12.5 


0.92 


11.0-14.1 


O.55M0/I, 


.385 


Hansen (Hurley) 
Hurley IFMR 


10-4/10-2 


Bergeron H 


12.4 


0.74 


> 10.5 


O.55M0/I, 


.385 


Hansen (default) 


10-4/10-2 


Hansen H 


12.1 


2.1 


10.9-12.5 


O.55M0/ • 




Hansen (default) 


10-4/10-2 
10-4/10-2 


Bergeron H 


12.4 


0.98 


10.5-14.1 


•••/••• 




Hansen (default) 


Bergeron H 


12.1 


0.76 


11.1-14.6 


•••/••• 










12.1 


0.74 


> 10.3 



Note. — The notation '• • •' indicates where we have marginalised over the parameter in question. 
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